1 /*
2 * read energy parameters from a file
3 *
4 * Stephan Kopp, Ivo Hofacker
5 * Vienna RNA Package
6 */
7
8 #ifdef HAVE_CONFIG_H
9 #include "config.h"
10 #endif
11
12 #include <stdio.h>
13 #include <stdlib.h>
14 #include <ctype.h>
15 #include <string.h>
16 #include <math.h>
17 #include <stdarg.h>
18 #include "ViennaRNA/utils/basic.h"
19 #include "ViennaRNA/io/utils.h"
20 #include "ViennaRNA/params/constants.h"
21 #include "ViennaRNA/params/default.h"
22 #include "ViennaRNA/params/io.h"
23 #include "ViennaRNA/static/energy_parameter_sets.h"
24
25
26 #define DEF -50
27 #define NST 0
28
29 /*
30 #################################
31 # PRIVATE VARIABLES #
32 #################################
33 */
34
35 PRIVATE char *last_param_file = NULL;
36
37 PRIVATE int stack_dim[2] = {
38 NBPAIRS + 1, NBPAIRS + 1
39 };
40 PRIVATE int stack_shift[2] = {
41 1, 1
42 };
43 PRIVATE int mismatch_dim[3] = {
44 NBPAIRS + 1, 5, 5
45 };
46 PRIVATE int mismatch_shift[3] = {
47 1, 0, 0
48 };
49 PRIVATE int int11_dim[4] = {
50 NBPAIRS + 1, NBPAIRS + 1, 5, 5
51 };
52 PRIVATE int int11_shift[4] = {
53 1, 1, 0, 0
54 };
55 PRIVATE int int21_dim[5] = {
56 NBPAIRS + 1, NBPAIRS + 1, 5, 5, 5
57 };
58 PRIVATE int int21_shift[5] = {
59 1, 1, 0, 0, 0
60 };
61 PRIVATE int int22_dim[6] = {
62 NBPAIRS + 1, NBPAIRS + 1, 5, 5, 5, 5
63 };
64 PRIVATE int int22_shift[6] = {
65 1, 1, 1, 1, 1, 1
66 };
67 PRIVATE int int22_post[6] = {
68 1, 1, 0, 0, 0, 0
69 };
70 PRIVATE int dangle_dim[2] = {
71 NBPAIRS + 1, 5
72 };
73 PRIVATE int dangle_shift[2] = {
74 1, 0
75 };
76
77 /*
78 #################################
79 # PRIVATE FUNCTION DECLARATIONS #
80 #################################
81 */
82
83 PRIVATE void
84 display_array(int *p,
85 int size,
86 int line,
87 FILE *fp);
88
89
90 PRIVATE char *
91 get_array1(char **content,
92 size_t *line_no,
93 int *arr,
94 int size);
95
96
97 PRIVATE void
98 ignore_comment(char *line);
99
100
101 PRIVATE void
102 check_symmetry(void);
103
104
105 PRIVATE void
106 update_nst(int array[NBPAIRS + 1][NBPAIRS + 1][5][5][5][5]);
107
108
109 /**
110 *** read a 1dimensional array from file
111 *** \param array a pointer to the first element in the array
112 *** \param dim the size of the array
113 *** \param shift the first position the new values will be written in
114 **/
115 PRIVATE void
116 rd_1dim(char **content,
117 size_t *line_no,
118 int *array,
119 int dim,
120 int shift);
121
122
123 PRIVATE void
124 rd_1dim_slice(char **content,
125 size_t *line_no,
126 int *array,
127 int dim,
128 int shift,
129 int post);
130
131
132 PRIVATE void
133 rd_2dim(char **content,
134 size_t *line_no,
135 int *array,
136 int dim[2],
137 int shift[2]);
138
139
140 PRIVATE void
141 rd_2dim_slice(char **content,
142 size_t *line_no,
143 int *array,
144 int dim[2],
145 int shift[2],
146 int post[2]);
147
148
149 PRIVATE void
150 rd_3dim(char **content,
151 size_t *line_no,
152 int *array,
153 int dim[3],
154 int shift[3]);
155
156
157 PRIVATE void
158 rd_3dim_slice(char **content,
159 size_t *line_no,
160 int *array,
161 int dim[3],
162 int shift[3],
163 int post[3]);
164
165
166 PRIVATE void
167 rd_4dim(char **content,
168 size_t *line_no,
169 int *array,
170 int dim[4],
171 int shift[4]);
172
173
174 PRIVATE void
175 rd_4dim_slice(char **content,
176 size_t *line_no,
177 int *array,
178 int dim[4],
179 int shift[4],
180 int post[4]);
181
182
183 PRIVATE void
184 rd_5dim(char **content,
185 size_t *line_no,
186 int *array,
187 int dim[5],
188 int shift[5]);
189
190
191 PRIVATE void
192 rd_5dim_slice(char **content,
193 size_t *line_no,
194 int *array,
195 int dim[5],
196 int shift[5],
197 int post[5]);
198
199
200 PRIVATE void
201 rd_6dim(char **content,
202 size_t *line_no,
203 int *array,
204 int dim[6],
205 int shift[6]);
206
207
208 PRIVATE void
209 rd_6dim_slice(char **content,
210 size_t *line_no,
211 int *array,
212 int dim[6],
213 int shift[6],
214 int post[6]);
215
216
217 PRIVATE void
218 rd_Tetraloop37(char **content,
219 size_t *line_no);
220
221
222 PRIVATE void
223 rd_Triloop37(char **content,
224 size_t *line_no);
225
226
227 PRIVATE void
228 rd_Hexaloop37(char **content,
229 size_t *line_no);
230
231
232 PRIVATE int
233 set_parameters_from_string(char **file_content,
234 const char *name);
235
236
237 PRIVATE char **
238 file2array(const char fname[]);
239
240
241 PRIVATE int
242 save_parameter_file(const char fname[],
243 unsigned int options);
244
245
246 /*
247 #################################
248 # BEGIN OF FUNCTION DEFINITIONS #
249 #################################
250 */
251 PUBLIC int
vrna_params_load(const char fname[],unsigned int options)252 vrna_params_load(const char fname[],
253 unsigned int options)
254 {
255 char *name, **file_content, **ptr;
256 int ret;
257
258 ret = 0;
259 file_content = file2array(fname);
260
261 if (file_content) {
262 name = vrna_basename(fname);
263
264 ret = set_parameters_from_string(file_content,
265 (const char *)name);
266
267 free(name);
268 for (ptr = file_content; *ptr != NULL; ptr++)
269 free(*ptr);
270
271 free(file_content);
272 }
273
274 return ret;
275 }
276
277
278 PUBLIC int
vrna_params_save(const char fname[],unsigned int options)279 vrna_params_save(const char fname[],
280 unsigned int options)
281 {
282 return save_parameter_file(fname, options);
283 }
284
285
286 PUBLIC int
vrna_params_load_from_string(const char * string,const char * name,unsigned int options)287 vrna_params_load_from_string(const char *string,
288 const char *name,
289 unsigned int options)
290 {
291 int ret = 0;
292
293 if (string) {
294 char **params_array, **ptr, *tmp_string, *token, *rest;
295 size_t lines, lines_mem;
296
297 lines = lines_mem = 0;
298 params_array = NULL;
299
300 /* convert string into array of lines */
301 tmp_string = strdup(string);
302
303 token = tmp_string;
304
305 /* tokenize the input string by separating lines at '\n' */
306 while (1) {
307 rest = strchr(token, '\n');
308 if (rest != NULL)
309 *rest = '\0';
310 else
311 break;
312
313 if (lines == lines_mem) {
314 lines_mem += 32768;
315 params_array = (char **)vrna_realloc(params_array, sizeof(char *) * lines_mem);
316 }
317
318 params_array[lines++] = strdup(token);
319
320 token = rest + 1;
321 }
322
323 /* reallocate to actual requirements */
324 params_array = (char **)vrna_realloc(params_array, sizeof(char *) * (lines + 1));
325 params_array[lines] = NULL;
326
327 /* actually apply parameters */
328 ret = set_parameters_from_string(params_array,
329 name);
330
331 /* cleanup memory */
332 free(tmp_string);
333
334 for (ptr = params_array; *ptr != NULL; ptr++)
335 free(*ptr);
336
337 free(params_array);
338 }
339
340 return ret;
341 }
342
343
344 PUBLIC int
vrna_params_load_defaults(void)345 vrna_params_load_defaults(void)
346 {
347 return vrna_params_load_RNA_Turner2004();
348 }
349
350
351 PUBLIC int
vrna_params_load_RNA_Turner2004(void)352 vrna_params_load_RNA_Turner2004(void)
353 {
354 return vrna_params_load_from_string((const char *)parameter_set_rna_turner2004,
355 "RNA - Turner 2004",
356 0);
357 }
358
359
360 PUBLIC int
vrna_params_load_RNA_Turner1999(void)361 vrna_params_load_RNA_Turner1999(void)
362 {
363 return vrna_params_load_from_string((const char *)parameter_set_rna_turner1999,
364 "RNA - Turner 1999",
365 0);
366 }
367
368
369 PUBLIC int
vrna_params_load_RNA_Andronescu2007(void)370 vrna_params_load_RNA_Andronescu2007(void)
371 {
372 return vrna_params_load_from_string((const char *)parameter_set_rna_andronescu2007,
373 "RNA - Andronescu 2007",
374 0);
375 }
376
377
378 PUBLIC int
vrna_params_load_RNA_Langdon2018(void)379 vrna_params_load_RNA_Langdon2018(void)
380 {
381 return vrna_params_load_from_string((const char *)parameter_set_rna_langdon2018,
382 "RNA - Langdon 2018",
383 0);
384 }
385
386
387 PUBLIC int
vrna_params_load_RNA_misc_special_hairpins(void)388 vrna_params_load_RNA_misc_special_hairpins(void)
389 {
390 return vrna_params_load_from_string((const char *)parameter_set_rna_misc_special_hairpins,
391 "RNA - Misc. Special Hairpins",
392 0);
393 }
394
395
396 PUBLIC int
vrna_params_load_DNA_Mathews2004(void)397 vrna_params_load_DNA_Mathews2004(void)
398 {
399 return vrna_params_load_from_string((const char *)parameter_set_dna_mathews2004,
400 "DNA - Mathews 2004",
401 0);
402 }
403
404
405 PUBLIC int
vrna_params_load_DNA_Mathews1999(void)406 vrna_params_load_DNA_Mathews1999(void)
407 {
408 return vrna_params_load_from_string((const char *)parameter_set_dna_mathews1999,
409 "DNA - Mathews 1999",
410 0);
411 }
412
413
414 /*
415 #####################################
416 # BEGIN OF STATIC HELPER FUNCTIONS #
417 #####################################
418 */
419 PRIVATE char **
file2array(const char fname[])420 file2array(const char fname[])
421 {
422 char **content, *line;
423 size_t lines_num, lines_mem;
424 FILE *fp;
425
426 content = NULL;
427
428 if ((fp = fopen(fname, "r"))) {
429 lines_num = 0;
430 lines_mem = 32768;
431
432 content = (char **)vrna_alloc(sizeof(char *) * lines_mem);
433
434 /* read file line-by-line */
435 while ((line = vrna_read_line(fp))) {
436 if (lines_num == lines_mem) {
437 lines_mem += 32768;
438 content = (char **)vrna_realloc(content, sizeof(char *) * lines_mem);
439 }
440
441 content[lines_num] = line;
442 lines_num++;
443 }
444
445 /* reallocate to actual requirements */
446 content = (char **)vrna_realloc(content, sizeof(char *) * (lines_num + 1));
447 content[lines_num] = NULL;
448
449 fclose(fp);
450 } else {
451 vrna_message_warning("read_parameter_file():"
452 "Can't open file %s\n",
453 fname);
454 }
455
456 return content;
457 }
458
459
460 PRIVATE int
set_parameters_from_string(char ** file_content,const char * name)461 set_parameters_from_string(char **file_content,
462 const char *name)
463 {
464 size_t line_no;
465 char *line, ident[256];
466 enum parset type;
467 int r;
468
469 line_no = 0;
470
471 if ((!file_content) ||
472 (!file_content[line_no]))
473 return 0;
474
475 /* store file name of parameter data set */
476 free(last_param_file);
477 last_param_file = (name) ? strdup(name) : NULL;
478
479 if (strncmp(file_content[line_no++], "## RNAfold parameter file v2.0", 30) != 0) {
480 vrna_message_warning("Missing header line in file.\n"
481 "May be this file has not v2.0 format.\n"
482 "Use INTERRUPT-key to stop.");
483 }
484
485 while ((line = file_content[line_no++])) {
486 r = sscanf(line, "# %255s", ident);
487 if (r == 1) {
488 type = gettype(ident);
489 switch (type) {
490 case QUIT:
491 break;
492 case S:
493 rd_2dim(file_content, &line_no, &(stack37[0][0]), stack_dim, stack_shift);
494 break;
495 case S_H:
496 rd_2dim(file_content, &line_no, &(stackdH[0][0]), stack_dim, stack_shift);
497 break;
498 case HP:
499 rd_1dim(file_content, &line_no, &(hairpin37[0]), 31, 0);
500 break;
501 case HP_H:
502 rd_1dim(file_content, &line_no, &(hairpindH[0]), 31, 0);
503 break;
504 case B:
505 rd_1dim(file_content, &line_no, &(bulge37[0]), 31, 0);
506 break;
507 case B_H:
508 rd_1dim(file_content, &line_no, &(bulgedH[0]), 31, 0);
509 break;
510 case IL:
511 rd_1dim(file_content, &line_no, &(internal_loop37[0]), 31, 0);
512 break;
513 case IL_H:
514 rd_1dim(file_content, &line_no, &(internal_loopdH[0]), 31, 0);
515 break;
516 case MME:
517 rd_3dim(file_content, &line_no, &(mismatchExt37[0][0][0]),
518 mismatch_dim,
519 mismatch_shift);
520 break;
521 case MME_H:
522 rd_3dim(file_content, &line_no, &(mismatchExtdH[0][0][0]),
523 mismatch_dim,
524 mismatch_shift);
525 break;
526 case MMH:
527 rd_3dim(file_content, &line_no, &(mismatchH37[0][0][0]),
528 mismatch_dim,
529 mismatch_shift);
530 break;
531 case MMH_H:
532 rd_3dim(file_content, &line_no, &(mismatchHdH[0][0][0]),
533 mismatch_dim,
534 mismatch_shift);
535 break;
536 case MMI:
537 rd_3dim(file_content, &line_no, &(mismatchI37[0][0][0]),
538 mismatch_dim,
539 mismatch_shift);
540 break;
541 case MMI_H:
542 rd_3dim(file_content, &line_no, &(mismatchIdH[0][0][0]),
543 mismatch_dim,
544 mismatch_shift);
545 break;
546 case MMI1N:
547 rd_3dim(file_content, &line_no, &(mismatch1nI37[0][0][0]),
548 mismatch_dim,
549 mismatch_shift);
550 break;
551 case MMI1N_H:
552 rd_3dim(file_content, &line_no, &(mismatch1nIdH[0][0][0]),
553 mismatch_dim,
554 mismatch_shift);
555 break;
556 case MMI23:
557 rd_3dim(file_content, &line_no, &(mismatch23I37[0][0][0]),
558 mismatch_dim,
559 mismatch_shift);
560 break;
561 case MMI23_H:
562 rd_3dim(file_content, &line_no, &(mismatch23IdH[0][0][0]),
563 mismatch_dim,
564 mismatch_shift);
565 break;
566 case MMM:
567 rd_3dim(file_content, &line_no, &(mismatchM37[0][0][0]),
568 mismatch_dim,
569 mismatch_shift);
570 break;
571 case MMM_H:
572 rd_3dim(file_content, &line_no, &(mismatchMdH[0][0][0]),
573 mismatch_dim,
574 mismatch_shift);
575 break;
576 case INT11:
577 rd_4dim(file_content, &line_no, &(int11_37[0][0][0][0]),
578 int11_dim,
579 int11_shift);
580 break;
581 case INT11_H:
582 rd_4dim(file_content, &line_no, &(int11_dH[0][0][0][0]),
583 int11_dim,
584 int11_shift);
585 break;
586 case INT21:
587 rd_5dim(file_content, &line_no, &(int21_37[0][0][0][0][0]),
588 int21_dim,
589 int21_shift);
590 break;
591 case INT21_H:
592 rd_5dim(file_content, &line_no, &(int21_dH[0][0][0][0][0]),
593 int21_dim,
594 int21_shift);
595 break;
596 case INT22:
597 rd_6dim_slice(file_content, &line_no, &(int22_37[0][0][0][0][0][0]),
598 int22_dim,
599 int22_shift,
600 int22_post);
601 update_nst(int22_37);
602 break;
603 case INT22_H:
604 rd_6dim_slice(file_content, &line_no, &(int22_dH[0][0][0][0][0][0]),
605 int22_dim,
606 int22_shift,
607 int22_post);
608 update_nst(int22_dH);
609 break;
610 case D5:
611 rd_2dim(file_content, &line_no, &(dangle5_37[0][0]),
612 dangle_dim,
613 dangle_shift);
614 break;
615 case D5_H:
616 rd_2dim(file_content, &line_no, &(dangle5_dH[0][0]),
617 dangle_dim,
618 dangle_shift);
619 break;
620 case D3:
621 rd_2dim(file_content, &line_no, &(dangle3_37[0][0]),
622 dangle_dim,
623 dangle_shift);
624 break;
625 case D3_H:
626 rd_2dim(file_content, &line_no, &(dangle3_dH[0][0]),
627 dangle_dim,
628 dangle_shift);
629 break;
630 case ML:
631 {
632 int values[6];
633 rd_1dim(file_content, &line_no, &values[0], 6, 0);
634 ML_BASE37 = values[0];
635 ML_BASEdH = values[1];
636 ML_closing37 = values[2];
637 ML_closingdH = values[3];
638 ML_intern37 = values[4];
639 ML_interndH = values[5];
640 }
641 break;
642 case NIN:
643 {
644 int values[3];
645 rd_1dim(file_content, &line_no, &values[0], 3, 0);
646 ninio37 = values[0];
647 niniodH = values[1];
648 MAX_NINIO = values[2];
649 }
650 break;
651 case MISC:
652 {
653 int values[4];
654 rd_1dim(file_content, &line_no, &values[0], 4, 0);
655 DuplexInit37 = values[0];
656 DuplexInitdH = values[1];
657 TerminalAU37 = values[2];
658 TerminalAUdH = values[3];
659 }
660 break;
661 case TL:
662 rd_Tetraloop37(file_content, &line_no);
663 break;
664 case TRI:
665 rd_Triloop37(file_content, &line_no);
666 break;
667 case HEX:
668 rd_Hexaloop37(file_content, &line_no);
669 break;
670 default: /* do nothing but complain */
671 vrna_message_warning("read_epars: Unknown field identifier in `%s'", line);
672 }
673 } /* else ignore line */
674 }
675
676 check_symmetry();
677 return 1;
678 }
679
680
681 /*------------------------------------------------------------*/
682
683 PRIVATE void
display_array(int * p,int size,int nl,FILE * fp)684 display_array(int *p,
685 int size,
686 int nl,
687 FILE *fp)
688 {
689 int i;
690
691 for (i = 1; i <= size; i++, p++) {
692 switch (*p) {
693 case INF:
694 fprintf(fp, " INF");
695 break;
696 case -INF:
697 fprintf(fp, " -INf");
698 break;
699 case DEF:
700 fprintf(fp, " DEF");
701 break;
702 default:
703 fprintf(fp, "%6d", *p);
704 break;
705 }
706 if ((i % nl) == 0)
707 fprintf(fp, "\n");
708 }
709 if (size % nl)
710 fprintf(fp, "\n");
711
712 return;
713 }
714
715
716 /*------------------------------------------------------------*/
717
718 PRIVATE char *
get_array1(char ** content,size_t * line_no,int * arr,int size)719 get_array1(char **content,
720 size_t *line_no,
721 int *arr,
722 int size)
723 {
724 int i, p, pos, pp, r, last;
725 char *line, buf[16];
726
727
728 i = last = 0;
729 while (i < size) {
730 line = content[(*line_no)++];
731 if (!line)
732 vrna_message_error("unexpected end of file in get_array1");
733
734 ignore_comment(line);
735 pos = 0;
736 while ((i < size) && (sscanf(line + pos, "%15s%n", buf, &pp) == 1)) {
737 pos += pp;
738 if (buf[0] == '*') {
739 i++;
740 continue;
741 } else if (buf[0] == 'x') {
742 /* should only be used for loop parameters */
743 if (i == 0)
744 vrna_message_error("can't extrapolate first value");
745
746 p = arr[last] + (int)(0.5 + lxc37 * log(((double)i) / (double)(last)));
747 } else if (strcmp(buf, "DEF") == 0) {
748 p = DEF;
749 } else if (strcmp(buf, "INF") == 0) {
750 p = INF;
751 } else if (strcmp(buf, "NST") == 0) {
752 p = NST;
753 } else {
754 r = sscanf(buf, "%d", &p);
755 if (r != 1) {
756 return line + pos;
757 vrna_message_error("can't interpret `%s' in get_array1", buf);
758 exit(1);
759 }
760
761 last = i;
762 }
763
764 arr[i++] = p;
765 }
766 }
767
768 return NULL;
769 }
770
771
772 PRIVATE void
rd_1dim(char ** content,size_t * line_no,int * array,int dim,int shift)773 rd_1dim(char **content,
774 size_t *line_no,
775 int *array,
776 int dim,
777 int shift)
778 {
779 rd_1dim_slice(content, line_no, array, dim, shift, 0);
780 }
781
782
783 PRIVATE void
rd_1dim_slice(char ** content,size_t * line_no,int * array,int dim,int shift,int post)784 rd_1dim_slice(char **content,
785 size_t *line_no,
786 int *array,
787 int dim,
788 int shift,
789 int post)
790 {
791 char *cp;
792
793 cp = get_array1(content, line_no, array + shift, dim - shift - post);
794
795 if (cp) {
796 vrna_message_error("\nrd_1dim: %s", cp);
797 exit(1);
798 }
799
800 return;
801 }
802
803
804 PRIVATE void
rd_2dim(char ** content,size_t * line_no,int * array,int dim[2],int shift[2])805 rd_2dim(char **content,
806 size_t *line_no,
807 int *array,
808 int dim[2],
809 int shift[2])
810 {
811 int post[2] = {
812 0, 0
813 };
814
815 rd_2dim_slice(content, line_no, array, dim, shift, post);
816 }
817
818
819 PRIVATE void
rd_2dim_slice(char ** content,size_t * line_no,int * array,int dim[2],int shift[2],int post[2])820 rd_2dim_slice(char **content,
821 size_t *line_no,
822 int *array,
823 int dim[2],
824 int shift[2],
825 int post[2])
826 {
827 int i;
828 int delta_pre = shift[0] + shift[1];
829 int delta_post = post[0] + post[1];
830
831 if (delta_pre + delta_post == 0) {
832 rd_1dim(content, line_no, array, dim[0] * dim[1], 0);
833 return;
834 }
835
836 for (i = shift[0]; i < dim[0] - post[0]; i++)
837 rd_1dim_slice(content, line_no, array + (i * dim[1]), dim[1], shift[1], post[1]);
838 return;
839 }
840
841
842 PRIVATE void
rd_3dim(char ** content,size_t * line_no,int * array,int dim[3],int shift[3])843 rd_3dim(char **content,
844 size_t *line_no,
845 int *array,
846 int dim[3],
847 int shift[3])
848 {
849 int post[3] = {
850 0, 0, 0
851 };
852
853 rd_3dim_slice(content, line_no, array,
854 dim,
855 shift,
856 post);
857 }
858
859
860 PRIVATE void
rd_3dim_slice(char ** content,size_t * line_no,int * array,int dim[3],int shift[3],int post[3])861 rd_3dim_slice(char **content,
862 size_t *line_no,
863 int *array,
864 int dim[3],
865 int shift[3],
866 int post[3])
867 {
868 int i;
869 int delta_pre = shift[0] + shift[1] + shift[2];
870 int delta_post = post[0] + post[1] + post[2];
871
872 if (delta_pre + delta_post == 0) {
873 rd_1dim(content, line_no, array, dim[0] * dim[1] * dim[2], 0);
874 return;
875 }
876
877 for (i = shift[0]; i < dim[0] - post[0]; i++) {
878 rd_2dim_slice(content, line_no, array + (i * dim[1] * dim[2]),
879 dim + 1,
880 shift + 1,
881 post + 1);
882 }
883 return;
884 }
885
886
887 PRIVATE void
rd_4dim(char ** content,size_t * line_no,int * array,int dim[4],int shift[4])888 rd_4dim(char **content,
889 size_t *line_no,
890 int *array,
891 int dim[4],
892 int shift[4])
893 {
894 int post[4] = {
895 0, 0, 0, 0
896 };
897
898 rd_4dim_slice(content, line_no, array,
899 dim,
900 shift,
901 post);
902 }
903
904
905 PRIVATE void
rd_4dim_slice(char ** content,size_t * line_no,int * array,int dim[4],int shift[4],int post[4])906 rd_4dim_slice(char **content,
907 size_t *line_no,
908 int *array,
909 int dim[4],
910 int shift[4],
911 int post[4])
912 {
913 int i;
914 int delta_pre = shift[0] + shift[1] + shift[2] + shift[3];
915 int delta_post = post[0] + post[1] + post[2] + post[3];
916
917 if (delta_pre + delta_post == 0) {
918 rd_1dim(content, line_no, array, dim[0] * dim[1] * dim[2] * dim[3], 0);
919 return;
920 }
921
922 for (i = shift[0]; i < dim[0] - post[0]; i++) {
923 rd_3dim_slice(content, line_no, array + (i * dim[1] * dim[2] * dim[3]),
924 dim + 1,
925 shift + 1,
926 post + 1);
927 }
928 return;
929 }
930
931
932 PRIVATE void
rd_5dim(char ** content,size_t * line_no,int * array,int dim[5],int shift[5])933 rd_5dim(char **content,
934 size_t *line_no,
935 int *array,
936 int dim[5],
937 int shift[5])
938 {
939 int post[5] = {
940 0, 0, 0, 0, 0
941 };
942
943 rd_5dim_slice(content, line_no, array,
944 dim,
945 shift,
946 post);
947 }
948
949
950 PRIVATE void
rd_5dim_slice(char ** content,size_t * line_no,int * array,int dim[5],int shift[5],int post[5])951 rd_5dim_slice(char **content,
952 size_t *line_no,
953 int *array,
954 int dim[5],
955 int shift[5],
956 int post[5])
957 {
958 int i;
959 int delta_pre = shift[0] + shift[1] + shift[2] + shift[3] + shift[4];
960 int delta_post = post[0] + post[1] + post[2] + post[3] + post[4];
961
962 if (delta_pre + delta_post == 0) {
963 rd_1dim(content, line_no, array, dim[0] * dim[1] * dim[2] * dim[3] * dim[4], 0);
964 return;
965 }
966
967 for (i = shift[0]; i < dim[0] - post[0]; i++)
968 rd_4dim_slice(content, line_no, array + (i * dim[1] * dim[2] * dim[3] * dim[4]),
969 dim + 1,
970 shift + 1,
971 post + 1);
972 return;
973 }
974
975
976 /**
977 *** \param dim1 The size of the first dimension
978 *** \param shift1 The pre shift for the first dimension
979 **/
980 PRIVATE void
rd_6dim(char ** content,size_t * line_no,int * array,int dim[6],int shift[6])981 rd_6dim(char **content,
982 size_t *line_no,
983 int *array,
984 int dim[6],
985 int shift[6])
986 {
987 int post[6] = {
988 0, 0, 0, 0, 0, 0
989 };
990
991 rd_6dim_slice(content, line_no, array,
992 dim,
993 shift,
994 post);
995 }
996
997
998 PRIVATE void
rd_6dim_slice(char ** content,size_t * line_no,int * array,int dim[6],int shift[6],int post[6])999 rd_6dim_slice(char **content,
1000 size_t *line_no,
1001 int *array,
1002 int dim[6],
1003 int shift[6],
1004 int post[6])
1005 {
1006 int i;
1007 int delta_pre = shift[0] + shift[1] + shift[2] + shift[3] + shift[4] + shift[5];
1008 int delta_post = post[0] + post[1] + post[2] + post[3] + post[4] + post[5];
1009
1010 if (delta_pre + delta_post == 0) {
1011 rd_1dim(content, line_no, array, dim[0] * dim[1] * dim[2] * dim[3] * dim[4] * dim[5], 0);
1012 return;
1013 }
1014
1015 for (i = shift[0]; i < dim[0] - post[0]; i++)
1016 rd_5dim_slice(content, line_no, array + (i * dim[1] * dim[2] * dim[3] * dim[4] * dim[5]),
1017 dim + 1,
1018 shift + 1,
1019 post + 1);
1020 return;
1021 }
1022
1023
1024 /*------------------------------------------------------------*/
1025 PRIVATE void
rd_Tetraloop37(char ** content,size_t * line_no)1026 rd_Tetraloop37(char **content,
1027 size_t *line_no)
1028 {
1029 int i, r;
1030 char *buf;
1031
1032 i = 0;
1033 /* erase old tetraloop entries */
1034 memset(&Tetraloops, 0, 281);
1035 memset(&Tetraloop37, 0, sizeof(int) * 40);
1036 memset(&TetraloopdH, 0, sizeof(int) * 40);
1037 do {
1038 buf = content[(*line_no)++];
1039 if (buf == NULL)
1040 break;
1041
1042 r = sscanf(buf, "%6s %d %d", &Tetraloops[7 * i], &Tetraloop37[i], &TetraloopdH[i]);
1043 strcat(Tetraloops, " ");
1044 i++;
1045 } while ((r == 3) && (i < 40));
1046
1047 /* decrease line number to continue parsing at line that failed before */
1048 (*line_no)--;
1049
1050 return;
1051 }
1052
1053
1054 /*------------------------------------------------------------*/
1055 PRIVATE void
rd_Hexaloop37(char ** content,size_t * line_no)1056 rd_Hexaloop37(char **content,
1057 size_t *line_no)
1058 {
1059 int i, r;
1060 char *buf;
1061
1062 i = 0;
1063 /* erase old hexaloop entries */
1064 memset(&Hexaloops, 0, 361);
1065 memset(&Hexaloop37, 0, sizeof(int) * 40);
1066 memset(&HexaloopdH, 0, sizeof(int) * 40);
1067 do {
1068 buf = content[(*line_no)++];
1069 if (buf == NULL)
1070 break;
1071
1072 r = sscanf(buf, "%8s %d %d", &Hexaloops[9 * i], &Hexaloop37[i], &HexaloopdH[i]);
1073 strcat(Hexaloops, " ");
1074 i++;
1075 } while ((r == 3) && (i < 40));
1076
1077 /* decrease line number to continue parsing at line that failed before */
1078 (*line_no)--;
1079
1080 return;
1081 }
1082
1083
1084 /*------------------------------------------------------------*/
1085 PRIVATE void
rd_Triloop37(char ** content,size_t * line_no)1086 rd_Triloop37(char **content,
1087 size_t *line_no)
1088 {
1089 int i, r;
1090 char *buf;
1091
1092 i = 0;
1093 /* erase old hexaloop entries */
1094 memset(&Triloops, 0, 241);
1095 memset(&Triloop37, 0, sizeof(int) * 40);
1096 memset(&TriloopdH, 0, sizeof(int) * 40);
1097 do {
1098 buf = content[(*line_no)++];
1099 if (buf == NULL)
1100 break;
1101
1102 r = sscanf(buf, "%5s %d %d", &Triloops[6 * i], &Triloop37[i], &TriloopdH[i]);
1103 strcat(Triloops, " ");
1104 i++;
1105 } while ((r == 3) && (i < 40));
1106
1107 /* decrease line number to continue parsing at line that failed before */
1108 (*line_no)--;
1109
1110 return;
1111 }
1112
1113
1114 /*------------------------------------------------------------*/
1115
1116
1117 PRIVATE void
ignore_comment(char * line)1118 ignore_comment(char *line)
1119 {
1120 /*
1121 * excise C style comments
1122 * only one comment per line, no multiline comments
1123 */
1124 char *cp1, *cp2;
1125
1126 if ((cp1 = strstr(line, "/*"))) {
1127 cp2 = strstr(cp1, "*/");
1128 if (cp2 == NULL)
1129 vrna_message_error("unclosed comment in parameter file");
1130
1131 /* can't use strcpy for overlapping strings */
1132 for (cp2 += 2; *cp2 != '\0'; cp2++, cp1++)
1133 *cp1 = *cp2;
1134 *cp1 = '\0';
1135 }
1136
1137 return;
1138 }
1139
1140
1141 /*------------------------------------------------------------*/
1142
1143 PRIVATE void
check_symmetry(void)1144 check_symmetry(void)
1145 {
1146 int i, j, k, l;
1147
1148 for (i = 0; i <= NBPAIRS; i++)
1149 for (j = 0; j <= NBPAIRS; j++)
1150 if (stack37[i][j] != stack37[j][i])
1151 vrna_message_warning("stacking energies not symmetric");
1152
1153 for (i = 0; i <= NBPAIRS; i++)
1154 for (j = 0; j <= NBPAIRS; j++)
1155 if (stackdH[i][j] != stackdH[j][i])
1156 vrna_message_warning("stacking enthalpies not symmetric");
1157
1158 /* interior 1x1 loops */
1159 for (i = 0; i <= NBPAIRS; i++)
1160 for (j = 0; j <= NBPAIRS; j++)
1161 for (k = 0; k < 5; k++)
1162 for (l = 0; l < 5; l++)
1163 if (int11_37[i][j][k][l] != int11_37[j][i][l][k])
1164 vrna_message_warning("int11 energies not symmetric (%d,%d,%d,%d) (%d vs. %d)",
1165 i, j, k, l, int11_37[i][j][k][l], int11_37[j][i][l][k]);
1166
1167 for (i = 0; i <= NBPAIRS; i++)
1168 for (j = 0; j <= NBPAIRS; j++)
1169 for (k = 0; k < 5; k++)
1170 for (l = 0; l < 5; l++)
1171 if (int11_dH[i][j][k][l] != int11_dH[j][i][l][k])
1172 vrna_message_warning("int11 enthalpies not symmetric");
1173
1174 /* interior 2x2 loops */
1175 for (i = 0; i <= NBPAIRS; i++)
1176 for (j = 0; j <= NBPAIRS; j++)
1177 for (k = 0; k < 5; k++)
1178 for (l = 0; l < 5; l++) {
1179 int m, n;
1180 for (m = 0; m < 5; m++)
1181 for (n = 0; n < 5; n++)
1182 if (int22_37[i][j][k][l][m][n] != int22_37[j][i][m][n][k][l])
1183 vrna_message_warning("int22 energies not symmetric");
1184 }
1185
1186 for (i = 0; i <= NBPAIRS; i++)
1187 for (j = 0; j <= NBPAIRS; j++)
1188 for (k = 0; k < 5; k++)
1189 for (l = 0; l < 5; l++) {
1190 int m, n;
1191 for (m = 0; m < 5; m++)
1192 for (n = 0; n < 5; n++)
1193 if (int22_dH[i][j][k][l][m][n] != int22_dH[j][i][m][n][k][l])
1194 vrna_message_warning("int22 enthalpies not symmetric: %d %d %d %d %d %d",
1195 i, j, k, l, m, n);
1196 }
1197 }
1198
1199
1200 /* update nonstandard nucleotide/basepair involved contributions for int22 */
1201 PRIVATE void
update_nst(int array[NBPAIRS+1][NBPAIRS+1][5][5][5][5])1202 update_nst(int array[NBPAIRS + 1][NBPAIRS + 1][5][5][5][5])
1203 {
1204 int i, j, k, l, m, n;
1205 int max, max2, max3, max4, max5, max6;
1206
1207 /* get maxima for one nonstandard nucleotide */
1208 for (i = 1; i < NBPAIRS; i++) {
1209 for (j = 1; j < NBPAIRS; j++) {
1210 for (k = 1; k < 5; k++) {
1211 for (l = 1; l < 5; l++) {
1212 for (m = 1; m < 5; m++) {
1213 max = max2 = max3 = max4 = -INF; /* max of {CGAU} */
1214 for (n = 1; n < 5; n++) {
1215 max = MAX2(max, array[i][j][k][l][m][n]);
1216 max2 = MAX2(max2, array[i][j][k][l][n][m]);
1217 max3 = MAX2(max3, array[i][j][k][n][l][m]);
1218 max4 = MAX2(max4, array[i][j][n][k][l][m]);
1219 }
1220 array[i][j][k][l][m][0] = max;
1221 array[i][j][k][l][0][m] = max2;
1222 array[i][j][k][0][l][m] = max3;
1223 array[i][j][0][k][l][m] = max4;
1224 }
1225 }
1226 }
1227 }
1228 }
1229 /* get maxima for two nonstandard nucleotides */
1230 for (i = 1; i < NBPAIRS; i++) {
1231 for (j = 1; j < NBPAIRS; j++) {
1232 for (k = 1; k < 5; k++) {
1233 for (l = 1; l < 5; l++) {
1234 max = max2 = max3 = max4 = max5 = max6 = -INF; /* max of {CGAU} */
1235 for (m = 1; m < 5; m++) {
1236 max = MAX2(max, array[i][j][k][l][m][0]);
1237 max2 = MAX2(max2, array[i][j][k][m][0][l]);
1238 max3 = MAX2(max3, array[i][j][m][0][k][l]);
1239 max4 = MAX2(max4, array[i][j][0][k][l][m]);
1240 max5 = MAX2(max5, array[i][j][0][k][m][l]);
1241 max6 = MAX2(max6, array[i][j][k][0][l][m]);
1242 }
1243 array[i][j][k][l][0][0] = max;
1244 array[i][j][k][0][0][l] = max2;
1245 array[i][j][0][0][k][l] = max3;
1246 array[i][j][k][0][l][0] = max6;
1247 array[i][j][0][k][0][l] = max5;
1248 array[i][j][0][k][l][0] = max4;
1249 }
1250 }
1251 }
1252 }
1253 /* get maxima for three nonstandard nucleotides */
1254 for (i = 1; i < NBPAIRS; i++) {
1255 for (j = 1; j < NBPAIRS; j++) {
1256 for (k = 1; k < 5; k++) {
1257 max = max2 = max3 = max4 = -INF; /* max of {CGAU} */
1258 for (l = 1; l < 5; l++) {
1259 /* should be arbitrary where index l resides in last 3 possible locations */
1260 max = MAX2(max, array[i][j][k][l][0][0]);
1261 max2 = MAX2(max2, array[i][j][0][k][l][0]);
1262 max3 = MAX2(max3, array[i][j][0][0][k][l]);
1263 max4 = MAX2(max4, array[i][j][0][0][l][k]);
1264 }
1265 array[i][j][k][0][0][0] = max;
1266 array[i][j][0][k][0][0] = max2;
1267 array[i][j][0][0][k][0] = max3;
1268 array[i][j][0][0][0][k] = max4;
1269 }
1270 }
1271 }
1272 /* get maxima for 4 nonstandard nucleotides */
1273 for (i = 1; i < NBPAIRS; i++) {
1274 for (j = 1; j < NBPAIRS; j++) {
1275 max = -INF; /* max of {CGAU} */
1276 for (k = 1; k < 5; k++)
1277 max = MAX2(max, array[i][j][k][0][0][0]);
1278 array[i][j][0][0][0][0] = max;
1279 }
1280 }
1281
1282 /*
1283 * now compute contributions for nonstandard base pairs ...
1284 * first, 1 nonstandard bp
1285 */
1286 for (i = 1; i < NBPAIRS; i++) {
1287 for (k = 0; k < 5; k++) {
1288 for (l = 0; l < 5; l++) {
1289 for (m = 0; m < 5; m++) {
1290 for (n = 0; n < 5; n++) {
1291 max = max2 = -INF;
1292 for (j = 1; j < NBPAIRS; j++) {
1293 max = MAX2(max, array[i][j][k][l][m][n]);
1294 max2 = MAX2(max2, array[j][i][k][l][m][n]);
1295 }
1296 array[i][NBPAIRS][k][l][m][n] = max;
1297 array[NBPAIRS][i][k][l][m][n] = max2;
1298 }
1299 }
1300 }
1301 }
1302 }
1303
1304 /* now 2 nst base pairs */
1305 for (k = 0; k < 5; k++) {
1306 for (l = 0; l < 5; l++) {
1307 for (m = 0; m < 5; m++) {
1308 for (n = 0; n < 5; n++) {
1309 max = -INF;
1310 for (j = 1; j < NBPAIRS; j++)
1311 max = MAX2(max, array[NBPAIRS][j][k][l][m][n]);
1312 array[NBPAIRS][NBPAIRS][k][l][m][n] = max;
1313 }
1314 }
1315 }
1316 }
1317 }
1318
1319
1320 PRIVATE int
save_parameter_file(const char fname[],unsigned int options)1321 save_parameter_file(const char fname[],
1322 unsigned int options)
1323 {
1324 FILE *outfp;
1325 int c;
1326 char *pnames[] = {
1327 "NP", "CG", "GC", "GU", "UG", "AU", "UA", " @"
1328 };
1329 char bnames[] = "@ACGU";
1330
1331 outfp = fopen(fname, "w");
1332 if (!outfp) {
1333 vrna_message_warning("can't open file %s", fname);
1334 return 0;
1335 }
1336
1337 fprintf(outfp, "## RNAfold parameter file v2.0\n");
1338
1339 fprintf(outfp, "\n# %s\n", settype(S));
1340 fprintf(outfp, "/* CG GC GU UG AU UA @ */\n");
1341 for (c = 1; c < NBPAIRS + 1; c++)
1342 display_array(stack37[c] + 1, NBPAIRS, NBPAIRS, outfp);
1343
1344 fprintf(outfp, "\n# %s\n", settype(S_H));
1345 fprintf(outfp, "/* CG GC GU UG AU UA @ */\n");
1346 for (c = 1; c < NBPAIRS + 1; c++)
1347 display_array(stackdH[c] + 1, NBPAIRS, NBPAIRS, outfp);
1348
1349 fprintf(outfp, "\n# %s\n", settype(MMH));
1350 {
1351 int i, k;
1352 for (k = 1; k < NBPAIRS + 1; k++)
1353 for (i = 0; i < 5; i++)
1354 display_array(mismatchH37[k][i], 5, 5, outfp);
1355 }
1356
1357 fprintf(outfp, "\n# %s\n", settype(MMH_H));
1358 {
1359 int i, k;
1360 for (k = 1; k < NBPAIRS + 1; k++)
1361 for (i = 0; i < 5; i++)
1362 display_array(mismatchHdH[k][i], 5, 5, outfp);
1363 }
1364
1365 fprintf(outfp, "\n# %s\n", settype(MMI));
1366 {
1367 int i, k;
1368 for (k = 1; k < NBPAIRS + 1; k++)
1369 for (i = 0; i < 5; i++)
1370 display_array(mismatchI37[k][i], 5, 5, outfp);
1371 }
1372
1373 fprintf(outfp, "\n# %s\n", settype(MMI_H));
1374 {
1375 int i, k;
1376 for (k = 1; k < NBPAIRS + 1; k++)
1377 for (i = 0; i < 5; i++)
1378 display_array(mismatchIdH[k][i], 5, 5, outfp);
1379 }
1380
1381 fprintf(outfp, "\n# %s\n", settype(MMI1N));
1382 {
1383 int i, k;
1384 for (k = 1; k < NBPAIRS + 1; k++)
1385 for (i = 0; i < 5; i++)
1386 display_array(mismatch1nI37[k][i], 5, 5, outfp);
1387 }
1388
1389 fprintf(outfp, "\n# %s\n", settype(MMI1N_H));
1390 {
1391 int i, k;
1392 for (k = 1; k < NBPAIRS + 1; k++)
1393 for (i = 0; i < 5; i++)
1394 display_array(mismatch1nIdH[k][i], 5, 5, outfp);
1395 }
1396
1397 fprintf(outfp, "\n# %s\n", settype(MMI23));
1398 {
1399 int i, k;
1400 for (k = 1; k < NBPAIRS + 1; k++)
1401 for (i = 0; i < 5; i++)
1402 display_array(mismatch23I37[k][i], 5, 5, outfp);
1403 }
1404
1405 fprintf(outfp, "\n# %s\n", settype(MMI23_H));
1406 {
1407 int i, k;
1408 for (k = 1; k < NBPAIRS + 1; k++)
1409 for (i = 0; i < 5; i++)
1410 display_array(mismatch23IdH[k][i], 5, 5, outfp);
1411 }
1412
1413 fprintf(outfp, "\n# %s\n", settype(MMM));
1414 {
1415 int i, k;
1416 for (k = 1; k < NBPAIRS + 1; k++)
1417 for (i = 0; i < 5; i++)
1418 display_array(mismatchM37[k][i], 5, 5, outfp);
1419 }
1420
1421 fprintf(outfp, "\n# %s\n", settype(MMM_H));
1422 {
1423 int i, k;
1424 for (k = 1; k < NBPAIRS + 1; k++)
1425 for (i = 0; i < 5; i++)
1426 display_array(mismatchMdH[k][i], 5, 5, outfp);
1427 }
1428
1429 fprintf(outfp, "\n# %s\n", settype(MME));
1430 {
1431 int i, k;
1432 for (k = 1; k < NBPAIRS + 1; k++)
1433 for (i = 0; i < 5; i++)
1434 display_array(mismatchExt37[k][i], 5, 5, outfp);
1435 }
1436
1437 fprintf(outfp, "\n# %s\n", settype(MME_H));
1438 {
1439 int i, k;
1440 for (k = 1; k < NBPAIRS + 1; k++)
1441 for (i = 0; i < 5; i++)
1442 display_array(mismatchExtdH[k][i], 5, 5, outfp);
1443 }
1444
1445 fprintf(outfp, "\n# %s\n", settype(D5));
1446 fprintf(outfp, "/* @ A C G U */\n");
1447 for (c = 1; c < NBPAIRS + 1; c++)
1448 display_array(dangle5_37[c], 5, 5, outfp);
1449
1450 fprintf(outfp, "\n# %s\n", settype(D5_H));
1451 fprintf(outfp, "/* @ A C G U */\n");
1452 for (c = 1; c < NBPAIRS + 1; c++)
1453 display_array(dangle5_dH[c], 5, 5, outfp);
1454
1455 fprintf(outfp, "\n# %s\n", settype(D3));
1456 fprintf(outfp, "/* @ A C G U */\n");
1457 for (c = 1; c < NBPAIRS + 1; c++)
1458 display_array(dangle3_37[c], 5, 5, outfp);
1459
1460 fprintf(outfp, "\n# %s\n", settype(D3_H));
1461 fprintf(outfp, "/* @ A C G U */\n");
1462 for (c = 1; c < NBPAIRS + 1; c++)
1463 display_array(dangle3_dH[c], 5, 5, outfp);
1464
1465
1466 /* dont print "no pair" entries for interior loop arrays */
1467 fprintf(outfp, "\n# %s\n", settype(INT11));
1468 {
1469 int i, k, l;
1470 for (k = 1; k < NBPAIRS + 1; k++)
1471 for (l = 1; l < NBPAIRS + 1; l++) {
1472 fprintf(outfp, "/* %2s..%2s */\n", pnames[k], pnames[l]);
1473 for (i = 0; i < 5; i++)
1474 display_array(int11_37[k][l][i], 5, 5, outfp);
1475 }
1476 }
1477
1478 fprintf(outfp, "\n# %s\n", settype(INT11_H));
1479 {
1480 int i, k, l;
1481 for (k = 1; k < NBPAIRS + 1; k++)
1482 for (l = 1; l < NBPAIRS + 1; l++) {
1483 fprintf(outfp, "/* %2s..%2s */\n", pnames[k], pnames[l]);
1484 for (i = 0; i < 5; i++)
1485 display_array(int11_dH[k][l][i], 5, 5, outfp);
1486 }
1487 }
1488
1489 fprintf(outfp, "\n# %s\n", settype(INT21));
1490 {
1491 int p1, p2, i, j;
1492 for (p1 = 1; p1 < NBPAIRS + 1; p1++)
1493 for (p2 = 1; p2 < NBPAIRS + 1; p2++)
1494 for (i = 0; i < 5; i++) {
1495 fprintf(outfp, "/* %2s.%c..%2s */\n",
1496 pnames[p1], bnames[i], pnames[p2]);
1497 for (j = 0; j < 5; j++)
1498 display_array(int21_37[p1][p2][i][j], 5, 5, outfp);
1499 }
1500 }
1501
1502 fprintf(outfp, "\n# %s\n", settype(INT21_H));
1503 {
1504 int p1, p2, i, j;
1505 for (p1 = 1; p1 < NBPAIRS + 1; p1++)
1506 for (p2 = 1; p2 < NBPAIRS + 1; p2++)
1507 for (i = 0; i < 5; i++) {
1508 fprintf(outfp, "/* %2s.%c..%2s */\n",
1509 pnames[p1], bnames[i], pnames[p2]);
1510 for (j = 0; j < 5; j++)
1511 display_array(int21_dH[p1][p2][i][j], 5, 5, outfp);
1512 }
1513 }
1514
1515 fprintf(outfp, "\n# %s\n", settype(INT22));
1516 {
1517 int p1, p2, i, j, k;
1518 for (p1 = 1; p1 < NBPAIRS; p1++)
1519 for (p2 = 1; p2 < NBPAIRS; p2++)
1520 for (i = 1; i < 5; i++)
1521 for (j = 1; j < 5; j++) {
1522 fprintf(outfp, "/* %2s.%c%c..%2s */\n",
1523 pnames[p1], bnames[i], bnames[j], pnames[p2]);
1524 for (k = 1; k < 5; k++)
1525 display_array(int22_37[p1][p2][i][j][k] + 1, 4, 5, outfp);
1526 }
1527 }
1528
1529 fprintf(outfp, "\n# %s\n", settype(INT22_H));
1530 {
1531 int p1, p2, i, j, k;
1532 for (p1 = 1; p1 < NBPAIRS; p1++)
1533 for (p2 = 1; p2 < NBPAIRS; p2++)
1534 for (i = 1; i < 5; i++)
1535 for (j = 1; j < 5; j++) {
1536 fprintf(outfp, "/* %2s.%c%c..%2s */\n",
1537 pnames[p1], bnames[i], bnames[j], pnames[p2]);
1538 for (k = 1; k < 5; k++)
1539 display_array(int22_dH[p1][p2][i][j][k] + 1, 4, 5, outfp);
1540 }
1541 }
1542
1543 fprintf(outfp, "\n# %s\n", settype(HP));
1544 display_array(hairpin37, 31, 10, outfp);
1545
1546 fprintf(outfp, "\n# %s\n", settype(HP_H));
1547 display_array(hairpindH, 31, 10, outfp);
1548
1549 fprintf(outfp, "\n# %s\n", settype(B));
1550 display_array(bulge37, 31, 10, outfp);
1551
1552 fprintf(outfp, "\n# %s\n", settype(B_H));
1553 display_array(bulgedH, 31, 10, outfp);
1554
1555 fprintf(outfp, "\n# %s\n", settype(IL));
1556 display_array(internal_loop37, 31, 10, outfp);
1557
1558 fprintf(outfp, "\n# %s\n", settype(IL_H));
1559 display_array(internal_loopdH, 31, 10, outfp);
1560
1561 fprintf(outfp, "\n# %s\n", settype(ML));
1562 fprintf(outfp, "/* F = cu*n_unpaired + cc + ci*loop_degree (+TermAU) */\n");
1563 fprintf(outfp, "/*\t cu\t cu_dH\t cc\t cc_dH\t ci\t ci_dH */\n");
1564 fprintf(outfp,
1565 "\t%6d\t%6d\t%6d\t%6d\t%6d\t%6d\n",
1566 ML_BASE37,
1567 ML_BASEdH,
1568 ML_closing37,
1569 ML_closingdH,
1570 ML_intern37,
1571 ML_interndH);
1572
1573 fprintf(outfp, "\n# %s\n", settype(NIN));
1574 fprintf(outfp, "/* Ninio = MIN(max, m*|n1-n2| */\n"
1575 "/*\t m\t m_dH max */\n"
1576 "\t%6d\t%6d\t%6d\n", ninio37, niniodH, MAX_NINIO);
1577
1578 fprintf(outfp, "\n# %s\n", settype(MISC));
1579 fprintf(outfp, "/* all parameters are pairs of 'energy enthalpy' */\n");
1580 fprintf(outfp, "/* DuplexInit TerminalAU LXC */\n");
1581 fprintf(outfp,
1582 " %6d %6d %6d %6d %3.6f %6d\n",
1583 DuplexInit37,
1584 DuplexInitdH,
1585 TerminalAU37,
1586 TerminalAUdH,
1587 lxc37,
1588 0);
1589
1590 fprintf(outfp, "\n# %s\n", settype(HEX));
1591 for (c = 0; c < strlen(Hexaloops) / 9; c++)
1592 fprintf(outfp, "\t%.8s %6d %6d\n", Hexaloops + c * 9, Hexaloop37[c], HexaloopdH[c]);
1593
1594 fprintf(outfp, "\n# %s\n", settype(TL));
1595 for (c = 0; c < strlen(Tetraloops) / 7; c++)
1596 fprintf(outfp, "\t%.6s %6d %6d\n", Tetraloops + c * 7, Tetraloop37[c], TetraloopdH[c]);
1597
1598 fprintf(outfp, "\n# %s\n", settype(TRI));
1599 for (c = 0; c < strlen(Triloops) / 6; c++)
1600 fprintf(outfp, "\t%.5s %6d %6d\n", Triloops + c * 6, Triloop37[c], TriloopdH[c]);
1601
1602 fprintf(outfp, "\n# %s\n", settype(QUIT));
1603 fclose(outfp);
1604
1605 return 1;
1606 }
1607
1608
1609 #ifndef VRNA_DISABLE_BACKWARD_COMPATIBILITY
1610
1611 /*
1612 *###########################################
1613 *# deprecated functions below #
1614 *###########################################
1615 */
1616 PUBLIC const char *
last_parameter_file(void)1617 last_parameter_file(void)
1618 {
1619 return (const char *)last_param_file;
1620 }
1621
1622
1623 /*------------------------------------------------------------*/
1624 PUBLIC void
read_parameter_file(const char fname[])1625 read_parameter_file(const char fname[])
1626 {
1627 (void)vrna_params_load(fname, VRNA_PARAMETER_FORMAT_DEFAULT);
1628 }
1629
1630
1631 PUBLIC char *
settype(enum parset s)1632 settype(enum parset s)
1633 {
1634 switch (s) {
1635 case S:
1636 return "stack";
1637 case S_H:
1638 return "stack_enthalpies";
1639 case HP:
1640 return "hairpin";
1641 case HP_H:
1642 return "hairpin_enthalpies";
1643 case B:
1644 return "bulge";
1645 case B_H:
1646 return "bulge_enthalpies";
1647 case IL:
1648 return "interior";
1649 case IL_H:
1650 return "interior_enthalpies";
1651 case MME:
1652 return "mismatch_exterior";
1653 case MME_H:
1654 return "mismatch_exterior_enthalpies";
1655 case MMH:
1656 return "mismatch_hairpin";
1657 case MMH_H:
1658 return "mismatch_hairpin_enthalpies";
1659 case MMI:
1660 return "mismatch_interior";
1661 case MMI_H:
1662 return "mismatch_interior_enthalpies";
1663 case MMI1N:
1664 return "mismatch_interior_1n";
1665 case MMI1N_H:
1666 return "mismatch_interior_1n_enthalpies";
1667 case MMI23:
1668 return "mismatch_interior_23";
1669 case MMI23_H:
1670 return "mismatch_interior_23_enthalpies";
1671 case MMM:
1672 return "mismatch_multi";
1673 case MMM_H:
1674 return "mismatch_multi_enthalpies";
1675 case D5:
1676 return "dangle5";
1677 case D5_H:
1678 return "dangle5_enthalpies";
1679 case D3:
1680 return "dangle3";
1681 case D3_H:
1682 return "dangle3_enthalpies";
1683 case INT11:
1684 return "int11";
1685 case INT11_H:
1686 return "int11_enthalpies";
1687 case INT21:
1688 return "int21";
1689 case INT21_H:
1690 return "int21_enthalpies";
1691 case INT22:
1692 return "int22";
1693 case INT22_H:
1694 return "int22_enthalpies";
1695 case ML:
1696 return "ML_params";
1697 case NIN:
1698 return "NINIO";
1699 case TRI:
1700 return "Triloops";
1701 case TL:
1702 return "Tetraloops";
1703 case HEX:
1704 return "Hexaloops";
1705 case QUIT:
1706 return "END";
1707 case MISC:
1708 return "Misc";
1709 default:
1710 vrna_message_error("\nThe answer is: 42\n");
1711 }
1712 return "";
1713 }
1714
1715
1716 /*------------------------------------------------------------*/
1717
1718 PUBLIC enum parset
gettype(const char * ident)1719 gettype(const char *ident)
1720 {
1721 if (strcmp(ident, "stack") == 0)
1722 return S;
1723 else if (strcmp(ident, "stack_enthalpies") == 0)
1724 return S_H;
1725 else if (strcmp(ident, "hairpin") == 0)
1726 return HP;
1727 else if (strcmp(ident, "hairpin_enthalpies") == 0)
1728 return HP_H;
1729 else if (strcmp(ident, "bulge") == 0)
1730 return B;
1731 else if (strcmp(ident, "bulge_enthalpies") == 0)
1732 return B_H;
1733 else if (strcmp(ident, "interior") == 0)
1734 return IL;
1735 else if (strcmp(ident, "interior_enthalpies") == 0)
1736 return IL_H;
1737 else if (strcmp(ident, "mismatch_exterior") == 0)
1738 return MME;
1739 else if (strcmp(ident, "mismatch_exterior_enthalpies") == 0)
1740 return MME_H;
1741 else if (strcmp(ident, "mismatch_hairpin") == 0)
1742 return MMH;
1743 else if (strcmp(ident, "mismatch_hairpin_enthalpies") == 0)
1744 return MMH_H;
1745 else if (strcmp(ident, "mismatch_interior") == 0)
1746 return MMI;
1747 else if (strcmp(ident, "mismatch_interior_enthalpies") == 0)
1748 return MMI_H;
1749 else if (strcmp(ident, "mismatch_interior_1n") == 0)
1750 return MMI1N;
1751 else if (strcmp(ident, "mismatch_interior_1n_enthalpies") == 0)
1752 return MMI1N_H;
1753 else if (strcmp(ident, "mismatch_interior_23") == 0)
1754 return MMI23;
1755 else if (strcmp(ident, "mismatch_interior_23_enthalpies") == 0)
1756 return MMI23_H;
1757 else if (strcmp(ident, "mismatch_multi") == 0)
1758 return MMM;
1759 else if (strcmp(ident, "mismatch_multi_enthalpies") == 0)
1760 return MMM_H;
1761 else if (strcmp(ident, "int11") == 0)
1762 return INT11;
1763 else if (strcmp(ident, "int11_enthalpies") == 0)
1764 return INT11_H;
1765 else if (strcmp(ident, "int21") == 0)
1766 return INT21;
1767 else if (strcmp(ident, "int21_enthalpies") == 0)
1768 return INT21_H;
1769 else if (strcmp(ident, "int22") == 0)
1770 return INT22;
1771 else if (strcmp(ident, "int22_enthalpies") == 0)
1772 return INT22_H;
1773 else if (strcmp(ident, "dangle5") == 0)
1774 return D5;
1775 else if (strcmp(ident, "dangle5_enthalpies") == 0)
1776 return D5_H;
1777 else if (strcmp(ident, "dangle3") == 0)
1778 return D3;
1779 else if (strcmp(ident, "dangle3_enthalpies") == 0)
1780 return D3_H;
1781 else if (strcmp(ident, "ML_params") == 0)
1782 return ML;
1783 else if (strcmp(ident, "NINIO") == 0)
1784 return NIN;
1785 else if (strcmp(ident, "Triloops") == 0)
1786 return TRI;
1787 else if (strcmp(ident, "Tetraloops") == 0)
1788 return TL;
1789 else if (strcmp(ident, "Hexaloops") == 0)
1790 return HEX;
1791 else if (strcmp(ident, "Misc") == 0)
1792 return MISC;
1793 else if (strcmp(ident, "END") == 0)
1794 return QUIT;
1795 else
1796 return UNKNOWN;
1797 }
1798
1799
1800 /*---------------------------------------------------------------*/
1801
1802 PUBLIC void
write_parameter_file(const char fname[])1803 write_parameter_file(const char fname[])
1804 {
1805 FILE *outfp;
1806 int c;
1807 char *pnames[] = {
1808 "NP", "CG", "GC", "GU", "UG", "AU", "UA", " @"
1809 };
1810 char bnames[] = "@ACGU";
1811
1812 outfp = fopen(fname, "w");
1813 if (!outfp) {
1814 vrna_message_error("can't open file %s", fname);
1815 exit(1);
1816 }
1817
1818 fprintf(outfp, "## RNAfold parameter file v2.0\n");
1819
1820 fprintf(outfp, "\n# %s\n", settype(S));
1821 fprintf(outfp, "/* CG GC GU UG AU UA @ */\n");
1822 for (c = 1; c < NBPAIRS + 1; c++)
1823 display_array(stack37[c] + 1, NBPAIRS, NBPAIRS, outfp);
1824
1825 fprintf(outfp, "\n# %s\n", settype(S_H));
1826 fprintf(outfp, "/* CG GC GU UG AU UA @ */\n");
1827 for (c = 1; c < NBPAIRS + 1; c++)
1828 display_array(stackdH[c] + 1, NBPAIRS, NBPAIRS, outfp);
1829
1830 fprintf(outfp, "\n# %s\n", settype(MMH));
1831 {
1832 int i, k;
1833 for (k = 1; k < NBPAIRS + 1; k++)
1834 for (i = 0; i < 5; i++)
1835 display_array(mismatchH37[k][i], 5, 5, outfp);
1836 }
1837
1838 fprintf(outfp, "\n# %s\n", settype(MMH_H));
1839 {
1840 int i, k;
1841 for (k = 1; k < NBPAIRS + 1; k++)
1842 for (i = 0; i < 5; i++)
1843 display_array(mismatchHdH[k][i], 5, 5, outfp);
1844 }
1845
1846 fprintf(outfp, "\n# %s\n", settype(MMI));
1847 {
1848 int i, k;
1849 for (k = 1; k < NBPAIRS + 1; k++)
1850 for (i = 0; i < 5; i++)
1851 display_array(mismatchI37[k][i], 5, 5, outfp);
1852 }
1853
1854 fprintf(outfp, "\n# %s\n", settype(MMI_H));
1855 {
1856 int i, k;
1857 for (k = 1; k < NBPAIRS + 1; k++)
1858 for (i = 0; i < 5; i++)
1859 display_array(mismatchIdH[k][i], 5, 5, outfp);
1860 }
1861
1862 fprintf(outfp, "\n# %s\n", settype(MMI1N));
1863 {
1864 int i, k;
1865 for (k = 1; k < NBPAIRS + 1; k++)
1866 for (i = 0; i < 5; i++)
1867 display_array(mismatch1nI37[k][i], 5, 5, outfp);
1868 }
1869
1870 fprintf(outfp, "\n# %s\n", settype(MMI1N_H));
1871 {
1872 int i, k;
1873 for (k = 1; k < NBPAIRS + 1; k++)
1874 for (i = 0; i < 5; i++)
1875 display_array(mismatch1nIdH[k][i], 5, 5, outfp);
1876 }
1877
1878 fprintf(outfp, "\n# %s\n", settype(MMI23));
1879 {
1880 int i, k;
1881 for (k = 1; k < NBPAIRS + 1; k++)
1882 for (i = 0; i < 5; i++)
1883 display_array(mismatch23I37[k][i], 5, 5, outfp);
1884 }
1885
1886 fprintf(outfp, "\n# %s\n", settype(MMI23_H));
1887 {
1888 int i, k;
1889 for (k = 1; k < NBPAIRS + 1; k++)
1890 for (i = 0; i < 5; i++)
1891 display_array(mismatch23IdH[k][i], 5, 5, outfp);
1892 }
1893
1894 fprintf(outfp, "\n# %s\n", settype(MMM));
1895 {
1896 int i, k;
1897 for (k = 1; k < NBPAIRS + 1; k++)
1898 for (i = 0; i < 5; i++)
1899 display_array(mismatchM37[k][i], 5, 5, outfp);
1900 }
1901
1902 fprintf(outfp, "\n# %s\n", settype(MMM_H));
1903 {
1904 int i, k;
1905 for (k = 1; k < NBPAIRS + 1; k++)
1906 for (i = 0; i < 5; i++)
1907 display_array(mismatchMdH[k][i], 5, 5, outfp);
1908 }
1909
1910 fprintf(outfp, "\n# %s\n", settype(MME));
1911 {
1912 int i, k;
1913 for (k = 1; k < NBPAIRS + 1; k++)
1914 for (i = 0; i < 5; i++)
1915 display_array(mismatchExt37[k][i], 5, 5, outfp);
1916 }
1917
1918 fprintf(outfp, "\n# %s\n", settype(MME_H));
1919 {
1920 int i, k;
1921 for (k = 1; k < NBPAIRS + 1; k++)
1922 for (i = 0; i < 5; i++)
1923 display_array(mismatchExtdH[k][i], 5, 5, outfp);
1924 }
1925
1926 fprintf(outfp, "\n# %s\n", settype(D5));
1927 fprintf(outfp, "/* @ A C G U */\n");
1928 for (c = 1; c < NBPAIRS + 1; c++)
1929 display_array(dangle5_37[c], 5, 5, outfp);
1930
1931 fprintf(outfp, "\n# %s\n", settype(D5_H));
1932 fprintf(outfp, "/* @ A C G U */\n");
1933 for (c = 1; c < NBPAIRS + 1; c++)
1934 display_array(dangle5_dH[c], 5, 5, outfp);
1935
1936 fprintf(outfp, "\n# %s\n", settype(D3));
1937 fprintf(outfp, "/* @ A C G U */\n");
1938 for (c = 1; c < NBPAIRS + 1; c++)
1939 display_array(dangle3_37[c], 5, 5, outfp);
1940
1941 fprintf(outfp, "\n# %s\n", settype(D3_H));
1942 fprintf(outfp, "/* @ A C G U */\n");
1943 for (c = 1; c < NBPAIRS + 1; c++)
1944 display_array(dangle3_dH[c], 5, 5, outfp);
1945
1946
1947 /* dont print "no pair" entries for interior loop arrays */
1948 fprintf(outfp, "\n# %s\n", settype(INT11));
1949 {
1950 int i, k, l;
1951 for (k = 1; k < NBPAIRS + 1; k++)
1952 for (l = 1; l < NBPAIRS + 1; l++) {
1953 fprintf(outfp, "/* %2s..%2s */\n", pnames[k], pnames[l]);
1954 for (i = 0; i < 5; i++)
1955 display_array(int11_37[k][l][i], 5, 5, outfp);
1956 }
1957 }
1958
1959 fprintf(outfp, "\n# %s\n", settype(INT11_H));
1960 {
1961 int i, k, l;
1962 for (k = 1; k < NBPAIRS + 1; k++)
1963 for (l = 1; l < NBPAIRS + 1; l++) {
1964 fprintf(outfp, "/* %2s..%2s */\n", pnames[k], pnames[l]);
1965 for (i = 0; i < 5; i++)
1966 display_array(int11_dH[k][l][i], 5, 5, outfp);
1967 }
1968 }
1969
1970 fprintf(outfp, "\n# %s\n", settype(INT21));
1971 {
1972 int p1, p2, i, j;
1973 for (p1 = 1; p1 < NBPAIRS + 1; p1++)
1974 for (p2 = 1; p2 < NBPAIRS + 1; p2++)
1975 for (i = 0; i < 5; i++) {
1976 fprintf(outfp, "/* %2s.%c..%2s */\n",
1977 pnames[p1], bnames[i], pnames[p2]);
1978 for (j = 0; j < 5; j++)
1979 display_array(int21_37[p1][p2][i][j], 5, 5, outfp);
1980 }
1981 }
1982
1983 fprintf(outfp, "\n# %s\n", settype(INT21_H));
1984 {
1985 int p1, p2, i, j;
1986 for (p1 = 1; p1 < NBPAIRS + 1; p1++)
1987 for (p2 = 1; p2 < NBPAIRS + 1; p2++)
1988 for (i = 0; i < 5; i++) {
1989 fprintf(outfp, "/* %2s.%c..%2s */\n",
1990 pnames[p1], bnames[i], pnames[p2]);
1991 for (j = 0; j < 5; j++)
1992 display_array(int21_dH[p1][p2][i][j], 5, 5, outfp);
1993 }
1994 }
1995
1996 fprintf(outfp, "\n# %s\n", settype(INT22));
1997 {
1998 int p1, p2, i, j, k;
1999 for (p1 = 1; p1 < NBPAIRS; p1++)
2000 for (p2 = 1; p2 < NBPAIRS; p2++)
2001 for (i = 1; i < 5; i++)
2002 for (j = 1; j < 5; j++) {
2003 fprintf(outfp, "/* %2s.%c%c..%2s */\n",
2004 pnames[p1], bnames[i], bnames[j], pnames[p2]);
2005 for (k = 1; k < 5; k++)
2006 display_array(int22_37[p1][p2][i][j][k] + 1, 4, 5, outfp);
2007 }
2008 }
2009
2010 fprintf(outfp, "\n# %s\n", settype(INT22_H));
2011 {
2012 int p1, p2, i, j, k;
2013 for (p1 = 1; p1 < NBPAIRS; p1++)
2014 for (p2 = 1; p2 < NBPAIRS; p2++)
2015 for (i = 1; i < 5; i++)
2016 for (j = 1; j < 5; j++) {
2017 fprintf(outfp, "/* %2s.%c%c..%2s */\n",
2018 pnames[p1], bnames[i], bnames[j], pnames[p2]);
2019 for (k = 1; k < 5; k++)
2020 display_array(int22_dH[p1][p2][i][j][k] + 1, 4, 5, outfp);
2021 }
2022 }
2023
2024 fprintf(outfp, "\n# %s\n", settype(HP));
2025 display_array(hairpin37, 31, 10, outfp);
2026
2027 fprintf(outfp, "\n# %s\n", settype(HP_H));
2028 display_array(hairpindH, 31, 10, outfp);
2029
2030 fprintf(outfp, "\n# %s\n", settype(B));
2031 display_array(bulge37, 31, 10, outfp);
2032
2033 fprintf(outfp, "\n# %s\n", settype(B_H));
2034 display_array(bulgedH, 31, 10, outfp);
2035
2036 fprintf(outfp, "\n# %s\n", settype(IL));
2037 display_array(internal_loop37, 31, 10, outfp);
2038
2039 fprintf(outfp, "\n# %s\n", settype(IL_H));
2040 display_array(internal_loopdH, 31, 10, outfp);
2041
2042 fprintf(outfp, "\n# %s\n", settype(ML));
2043 fprintf(outfp, "/* F = cu*n_unpaired + cc + ci*loop_degree (+TermAU) */\n");
2044 fprintf(outfp, "/*\t cu\t cu_dH\t cc\t cc_dH\t ci\t ci_dH */\n");
2045 fprintf(outfp,
2046 "\t%6d\t%6d\t%6d\t%6d\t%6d\t%6d\n",
2047 ML_BASE37,
2048 ML_BASEdH,
2049 ML_closing37,
2050 ML_closingdH,
2051 ML_intern37,
2052 ML_interndH);
2053
2054 fprintf(outfp, "\n# %s\n", settype(NIN));
2055 fprintf(outfp, "/* Ninio = MIN(max, m*|n1-n2| */\n"
2056 "/*\t m\t m_dH max */\n"
2057 "\t%6d\t%6d\t%6d\n", ninio37, niniodH, MAX_NINIO);
2058
2059 fprintf(outfp, "\n# %s\n", settype(MISC));
2060 fprintf(outfp, "/* all parameters are pairs of 'energy enthalpy' */\n");
2061 fprintf(outfp, "/* DuplexInit TerminalAU LXC */\n");
2062 fprintf(outfp,
2063 " %6d %6d %6d %6d %3.6f %6d\n",
2064 DuplexInit37,
2065 DuplexInitdH,
2066 TerminalAU37,
2067 TerminalAUdH,
2068 lxc37,
2069 0);
2070
2071 fprintf(outfp, "\n# %s\n", settype(HEX));
2072 for (c = 0; c < strlen(Hexaloops) / 9; c++)
2073 fprintf(outfp, "\t%.8s %6d %6d\n", Hexaloops + c * 9, Hexaloop37[c], HexaloopdH[c]);
2074
2075 fprintf(outfp, "\n# %s\n", settype(TL));
2076 for (c = 0; c < strlen(Tetraloops) / 7; c++)
2077 fprintf(outfp, "\t%.6s %6d %6d\n", Tetraloops + c * 7, Tetraloop37[c], TetraloopdH[c]);
2078
2079 fprintf(outfp, "\n# %s\n", settype(TRI));
2080 for (c = 0; c < strlen(Triloops) / 6; c++)
2081 fprintf(outfp, "\t%.5s %6d %6d\n", Triloops + c * 6, Triloop37[c], TriloopdH[c]);
2082
2083 fprintf(outfp, "\n# %s\n", settype(QUIT));
2084 fclose(outfp);
2085 }
2086
2087
2088 #endif
2089