--- /dev/null
+/*
+ read energy parameters from a file
+
+ Stephan Kopp, Ivo Hofacker
+ Vienna RNA Package
+*/
+
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <ctype.h>
+#include <string.h>
+#include <math.h>
+#include <stdarg.h>
+#include "utils.h"
+#include "energy_const.h"
+#include "energy_par.h"
+#include "read_epars.h"
+
+#define PUBLIC
+#define PRIVATE static
+#define PARSET 20
+
+#define DEF -50
+#define NST 0
+
+PRIVATE FILE *fp;
+
+PRIVATE void display_array(int *p, int size, int line, FILE *fp);
+PRIVATE char *get_array1(int *arr, int size);
+PRIVATE void ignore_comment(char *line);
+PRIVATE void check_symmetry(void);
+PRIVATE void update_nst(int array[NBPAIRS+1][NBPAIRS+1][5][5][5][5]);
+
+/**
+*** read a 1dimensional array from file
+*** \param array a pointer to the first element in the array
+*** \param dim the size of the array
+*** \param shift the first position the new values will be written in
+**/
+PRIVATE void rd_1dim(int *array, int dim, int shift);
+PRIVATE void rd_1dim_slice(int *array, int dim, int shift, int post);
+PRIVATE void rd_2dim(int *array,
+ int dim1, int dim2,
+ int shift1, int shift2);
+PRIVATE void rd_2dim_slice(int *array,
+ int dim1, int dim2,
+ int shift1, int shift2,
+ int post1, int post2);
+PRIVATE void rd_3dim(int *array,
+ int dim1, int dim2, int dim3,
+ int shift1, int shift2, int shift3);
+PRIVATE void rd_3dim_slice(int *array,
+ int dim1, int dim2, int dim3,
+ int shift1, int shift2, int shift3,
+ int post1, int post2, int post3);
+PRIVATE void rd_4dim(int *array,
+ int dim1, int dim2, int dim3, int dim4,
+ int shift1, int shift2, int shift3, int shift4);
+PRIVATE void rd_4dim_slice(int *array,
+ int dim1, int dim2, int dim3, int dim4,
+ int shift1, int shift2, int shift3, int shift4,
+ int post1, int post2, int post3, int post4);
+PRIVATE void rd_5dim(int *array,
+ int dim1, int dim2, int dim3, int dim4, int dim5,
+ int shift1, int shift2, int shift3, int shift4, int shift5);
+PRIVATE void rd_5dim_slice(int *array,
+ int dim1, int dim2, int dim3, int dim4, int dim5,
+ int shift1, int shift2, int shift3, int shift4, int shift5,
+ int post1, int post2, int post3, int post4, int post5);
+PRIVATE void rd_6dim(int *array,
+ int dim1, int dim2, int dim3, int dim4, int dim5, int dim6,
+ int shift1, int shift2, int shift3, int shift4, int shift5, int shift6);
+PRIVATE void rd_6dim_slice(int *array,
+ int dim1, int dim2, int dim3, int dim4, int dim5, int dim6,
+ int shift1, int shift2, int shift3, int shift4, int shift5, int shift6,
+ int post1, int post2, int post3, int post4, int post5, int post6);
+PRIVATE void rd_Tetraloop37(void);
+PRIVATE void rd_Triloop37(void);
+PRIVATE void rd_Hexaloop37(void);
+
+/*------------------------------------------------------------*/
+PUBLIC void read_parameter_file(const char fname[]){
+ char *line, ident[256];
+ enum parset type;
+ int r;
+
+ if (!(fp=fopen(fname,"r"))) {
+ fprintf(stderr,
+ "\nread_parameter_file:\n"
+ "\t\tcan't open file %s\n"
+ "\t\tusing default parameters instead.\n", fname);
+ return;
+ }
+
+ if (!(line = get_line(fp))) {
+ fprintf(stderr," File %s is inproper.\n", fname);
+ fclose(fp);
+ return;
+ }
+
+ if (strncmp(line,"## RNAfold parameter file v2.0",30)!=0) {
+ fprintf(stderr,
+ "Missing header line in file.\n"
+ "May be this file has not v2.0 format.\n"
+ "Use INTERRUPT-key to stop.\n");
+ }
+ free(line);
+
+ while((line=get_line(fp))) {
+
+ r = sscanf(line, "# %255s", ident);
+ if (r==1) {
+ type = gettype(ident);
+ switch (type){
+ case QUIT: break;
+ case S: rd_2dim(&(stack37[0][0]), NBPAIRS+1, NBPAIRS+1, 1, 1);
+ break;
+ case S_H: rd_2dim(&(stackdH[0][0]), NBPAIRS+1, NBPAIRS+1, 1, 1);
+ break;
+ case HP: rd_1dim(&(hairpin37[0]), 31, 0);
+ break;
+ case HP_H: rd_1dim(&(hairpindH[0]), 31, 0);
+ break;
+ case B: rd_1dim(&(bulge37[0]), 31, 0);
+ break;
+ case B_H: rd_1dim(&(bulgedH[0]), 31, 0);
+ break;
+ case IL: rd_1dim(&(internal_loop37[0]), 31, 0);
+ break;
+ case IL_H: rd_1dim(&(internal_loopdH[0]), 31, 0);
+ break;
+ case MME: rd_3dim(&(mismatchExt37[0][0][0]),
+ NBPAIRS+1, 5, 5,
+ 1, 0, 0);
+ break;
+ case MME_H: rd_3dim(&(mismatchExtdH[0][0][0]),
+ NBPAIRS+1, 5, 5,
+ 1, 0, 0);
+ break;
+ case MMH: rd_3dim(&(mismatchH37[0][0][0]),
+ NBPAIRS+1, 5, 5,
+ 1, 0, 0);
+ break;
+ case MMH_H: rd_3dim(&(mismatchHdH[0][0][0]),
+ NBPAIRS+1, 5, 5,
+ 1, 0, 0);
+ break;
+ case MMI: rd_3dim(&(mismatchI37[0][0][0]),
+ NBPAIRS+1, 5, 5,
+ 1, 0, 0);
+ break;
+ case MMI_H: rd_3dim(&(mismatchIdH[0][0][0]),
+ NBPAIRS+1, 5, 5,
+ 1, 0, 0);
+ break;
+ case MMI1N: rd_3dim(&(mismatch1nI37[0][0][0]),
+ NBPAIRS+1, 5, 5,
+ 1, 0, 0);
+ break;
+ case MMI1N_H: rd_3dim(&(mismatch1nIdH[0][0][0]),
+ NBPAIRS+1, 5, 5,
+ 1, 0, 0);
+ break;
+ case MMI23: rd_3dim(&(mismatch23I37[0][0][0]),
+ NBPAIRS+1, 5, 5,
+ 1, 0, 0);
+ break;
+ case MMI23_H: rd_3dim(&(mismatch23IdH[0][0][0]),
+ NBPAIRS+1, 5, 5,
+ 1, 0, 0);
+ break;
+ case MMM: rd_3dim(&(mismatchM37[0][0][0]),
+ NBPAIRS+1, 5, 5,
+ 1, 0, 0);
+ break;
+ case MMM_H: rd_3dim(&(mismatchMdH[0][0][0]),
+ NBPAIRS+1, 5, 5,
+ 1, 0, 0);
+ break;
+ case INT11: rd_4dim(&(int11_37[0][0][0][0]),
+ NBPAIRS+1, NBPAIRS+1, 5, 5,
+ 1, 1, 0, 0);
+ break;
+ case INT11_H: rd_4dim(&(int11_dH[0][0][0][0]),
+ NBPAIRS+1, NBPAIRS+1, 5, 5,
+ 1, 1, 0, 0);
+ break;
+ case INT21: rd_5dim(&(int21_37[0][0][0][0][0]),
+ NBPAIRS+1, NBPAIRS+1, 5, 5, 5,
+ 1, 1, 0, 0, 0);
+ break;
+ case INT21_H: rd_5dim(&(int21_dH[0][0][0][0][0]),
+ NBPAIRS+1, NBPAIRS+1, 5, 5, 5,
+ 1, 1, 0, 0, 0);
+ break;
+ case INT22: rd_6dim_slice(&(int22_37[0][0][0][0][0][0]),
+ NBPAIRS+1, NBPAIRS+1, 5, 5, 5, 5,
+ 1, 1, 1, 1, 1, 1,
+ 1, 1, 0, 0, 0, 0);
+ update_nst(int22_37);
+ break;
+ case INT22_H: rd_6dim_slice(&(int22_dH[0][0][0][0][0][0]),
+ NBPAIRS+1, NBPAIRS+1, 5, 5, 5, 5,
+ 1, 1, 1, 1, 1, 1,
+ 1, 1, 0, 0, 0, 0);
+ update_nst(int22_dH);
+ break;
+ case D5: rd_2dim(&(dangle5_37[0][0]), NBPAIRS+1, 5, 1, 0);
+ break;
+ case D5_H: rd_2dim(&(dangle5_dH[0][0]), NBPAIRS+1, 5, 1, 0);
+ break;
+ case D3: rd_2dim(&(dangle3_37[0][0]), NBPAIRS+1, 5, 1, 0);
+ break;
+ case D3_H: rd_2dim(&(dangle3_dH[0][0]), NBPAIRS+1, 5, 1, 0);
+ break;
+ case ML: {
+ int values[6];
+ rd_1dim(&values[0], 6, 0);
+ ML_BASE37 = values[0];
+ ML_BASEdH = values[1];
+ ML_closing37 = values[2];
+ ML_closingdH = values[3];
+ ML_intern37 = values[4];
+ ML_interndH = values[5];
+ }
+ break;
+ case NIN: {
+ int values[3];
+ rd_1dim(&values[0], 3, 0);
+ ninio37 = values[0];
+ niniodH = values[1];
+ MAX_NINIO = values[2];
+ }
+ break;
+ case MISC: {
+ int values[4];
+ rd_1dim(&values[0], 4, 0);
+ DuplexInit37 = values[0];
+ DuplexInitdH = values[1];
+ TerminalAU37 = values[2];
+ TerminalAUdH = values[3];
+ }
+ break;
+ case TL: rd_Tetraloop37();
+ break;
+ case TRI: rd_Triloop37();
+ break;
+ case HEX: rd_Hexaloop37();
+ break;
+ default: /* do nothing but complain */
+ fprintf(stderr,"read_epars: Unknown field identifier in `%s'\n", line);
+ }
+ } /* else ignore line */
+ free(line);
+ }
+ fclose(fp);
+
+ check_symmetry();
+ return;
+}
+
+/*------------------------------------------------------------*/
+
+PRIVATE void display_array(int *p, int size, int nl, FILE *fp)
+{
+ int i;
+
+ for (i=1; i<=size; i++, p++) {
+ switch(*p)
+ {
+ case INF: fprintf(fp," INF"); break;
+ case -INF: fprintf(fp," -INf"); break;
+ case DEF: fprintf(fp," DEF"); break;
+ default: fprintf(fp,"%6d", *p); break;
+ }
+ if ((i%nl)==0) fprintf(fp,"\n");
+ }
+ if (size%nl) fprintf(fp,"\n");
+ return;
+}
+
+/*------------------------------------------------------------*/
+
+PRIVATE char *get_array1(int *arr, int size)
+{
+ int i, p, pos, pp, r, last;
+ char *line, buf[16];
+
+
+ i = last = 0;
+ while( i<size ) {
+ line = get_line(fp);
+ if (!line) nrerror("unexpected end of file in get_array1");
+ ignore_comment(line);
+ pos=0;
+ while ((i<size)&&(sscanf(line+pos,"%15s%n", buf, &pp)==1)) {
+ pos += pp;
+ if (buf[0]=='*') {i++; continue;}
+ else if (buf[0]=='x') { /* should only be used for loop parameters */
+ if (i==0) nrerror("can't extrapolate first value");
+ p = arr[last] + (int) (0.5+ lxc37*log(((double) i)/(double)(last)));
+ }
+ else if (strcmp(buf,"DEF") == 0) p = DEF;
+ else if (strcmp(buf,"INF") == 0) p = INF;
+ else if (strcmp(buf,"NST") == 0) p = NST;
+ else {
+ r=sscanf(buf,"%d", &p);
+ if (r!=1) {
+ return line+pos;
+ fprintf(stderr, "can't interpret `%s' in get_array1\n", buf);
+ exit(1);
+ }
+ last = i;
+ }
+ arr[i++]=p;
+ }
+ free(line);
+ }
+
+ return NULL;
+}
+
+PRIVATE void rd_1dim(int *array, int dim, int shift){
+ rd_1dim_slice(array, dim, shift, 0);
+}
+
+PRIVATE void rd_1dim_slice(int *array, int dim, int shift, int post){
+ char *cp;
+ cp = get_array1(array+shift, dim-shift-post);
+
+ if (cp) {
+ fprintf(stderr,"\nrd_1dim: %s\n", cp);
+ exit(1);
+ }
+ return;
+}
+
+PRIVATE void rd_2dim(int *array, int dim1, int dim2, int shift1, int shift2){
+ rd_2dim_slice(array, dim1, dim2, shift1, shift2, 0, 0);
+}
+
+PRIVATE void rd_2dim_slice(int *array,
+ int dim1, int dim2,
+ int shift1, int shift2,
+ int post1, int post2){
+ int i;
+ int delta_pre = shift1 + shift2;
+ int delta_post = post1 + post2;
+
+ if(delta_pre + delta_post == 0){
+ rd_1dim(array, dim1 * dim2, 0);
+ return;
+ }
+ for (i=shift1; i<dim1 - post1; i++)
+ rd_1dim_slice(array + (i*dim2), dim2, shift2, post2);
+ return;
+}
+
+PRIVATE void rd_3dim(int *array, int dim1, int dim2, int dim3, int shift1, int shift2, int shift3){
+ rd_3dim_slice(array,
+ dim1, dim2, dim3,
+ shift1, shift2, shift3,
+ 0, 0, 0);
+}
+
+PRIVATE void rd_3dim_slice(int *array,
+ int dim1, int dim2, int dim3,
+ int shift1, int shift2, int shift3,
+ int post1, int post2, int post3){
+ int i;
+ int delta_pre = shift1 + shift2 + shift3;
+ int delta_post = post1 + post2 + post3;
+
+ if(delta_pre + delta_post == 0){
+ rd_1dim(array, dim1 * dim2 * dim3, 0);
+ return;
+ }
+ for (i=shift1; i<dim1 - post1; i++){
+ rd_2dim_slice(array + (i * dim2 * dim3),
+ dim2, dim3,
+ shift2, shift3,
+ post2, post3);
+ }
+ return;
+}
+
+PRIVATE void rd_4dim(int *array,
+ int dim1, int dim2, int dim3, int dim4,
+ int shift1, int shift2, int shift3, int shift4){
+ rd_4dim_slice(array,
+ dim1, dim2, dim3, dim4,
+ shift1, shift2, shift3, shift4,
+ 0, 0, 0, 0);
+}
+
+PRIVATE void rd_4dim_slice(int *array,
+ int dim1, int dim2, int dim3, int dim4,
+ int shift1, int shift2, int shift3, int shift4,
+ int post1, int post2, int post3, int post4){
+ int i;
+ int delta_pre = shift1 + shift2 + shift3 + shift4;
+ int delta_post = post1 + post2 + post3 + post4;
+
+ if(delta_pre + delta_post == 0){
+ rd_1dim(array, dim1 * dim2 * dim3 * dim4, 0);
+ return;
+ }
+ for(i=shift1; i<dim1 - post1; i++){
+ rd_3dim_slice(array + (i * dim2 * dim3 * dim4),
+ dim2, dim3, dim4,
+ shift2, shift3, shift4,
+ post2, post3, post4);
+ }
+ return;
+}
+
+PRIVATE void rd_5dim(int *array,
+ int dim1, int dim2, int dim3, int dim4, int dim5,
+ int shift1, int shift2, int shift3, int shift4, int shift5){
+ rd_5dim_slice(array,
+ dim1, dim2, dim3, dim4, dim5,
+ shift1, shift2, shift3, shift4, shift5,
+ 0, 0, 0, 0, 0);
+}
+
+PRIVATE void rd_5dim_slice(int *array,
+ int dim1, int dim2, int dim3, int dim4, int dim5,
+ int shift1, int shift2, int shift3, int shift4, int shift5,
+ int post1, int post2, int post3, int post4, int post5){
+ int i;
+ int delta_pre = shift1 + shift2 + shift3 + shift4 + shift5;
+ int delta_post = post1 + post2 + post3 + post4 + post5;
+
+ if(delta_pre + delta_post == 0){
+ rd_1dim(array, dim1 * dim2 * dim3 * dim4 * dim5, 0);
+ return;
+ }
+ for(i=shift1; i<dim1 - post1; i++)
+ rd_4dim_slice(array + (i * dim2 * dim3 * dim4 * dim5),
+ dim2, dim3, dim4, dim5,
+ shift2, shift3, shift4, shift5,
+ post2, post3, post4, post5);
+ return;
+}
+
+/**
+*** \param dim1 The size of the first dimension
+*** \param shift1 The pre shift for the first dimension
+**/
+PRIVATE void rd_6dim(int *array,
+ int dim1, int dim2, int dim3, int dim4, int dim5, int dim6,
+ int shift1, int shift2, int shift3, int shift4, int shift5, int shift6){
+ rd_6dim_slice(array,
+ dim1, dim2, dim3, dim4, dim5, dim6,
+ shift1, shift2, shift3, shift4, shift5, shift6,
+ 0, 0, 0, 0, 0, 0);
+}
+
+PRIVATE void rd_6dim_slice(int *array,
+ int dim1, int dim2, int dim3, int dim4, int dim5, int dim6,
+ int shift1, int shift2, int shift3, int shift4, int shift5, int shift6,
+ int post1, int post2, int post3, int post4, int post5, int post6){
+ int i;
+ int delta_pre = shift1 + shift2 + shift3 + shift4 + shift5 + shift6;
+ int delta_post = post1 + post2 + post3 + post4 + post5 + post6;
+
+ if(delta_pre + delta_post == 0){
+ rd_1dim(array, dim1 * dim2 * dim3 * dim4 * dim5 * dim6, 0);
+ return;
+ }
+ for(i=shift1; i<dim1 - post1; i++)
+ rd_5dim_slice(array + (i * dim2 * dim3 * dim4 * dim5 * dim6),
+ dim2, dim3, dim4, dim5, dim6,
+ shift2, shift3, shift4, shift5, shift6,
+ post2, post3, post4, post5, post6);
+ return;
+}
+
+
+/*------------------------------------------------------------*/
+PRIVATE void rd_Tetraloop37(void)
+{
+ int i, r;
+ char *buf;
+
+ i=0;
+ /* erase old tetraloop entries */
+ memset(&Tetraloops, 0, 281);
+ memset(&Tetraloop37, 0, sizeof(int)*40);
+ memset(&TetraloopdH, 0, sizeof(int)*40);
+ do {
+ buf = get_line(fp);
+ if (buf==NULL) break;
+ r = sscanf(buf,"%6s %d %d", &Tetraloops[7*i], &Tetraloop37[i], &TetraloopdH[i]);
+ strcat(Tetraloops, " ");
+ free(buf);
+ i++;
+ } while((r==3)&&(i<40));
+ return;
+}
+
+/*------------------------------------------------------------*/
+PRIVATE void rd_Hexaloop37(void)
+{
+ int i, r;
+ char *buf;
+
+ i=0;
+ /* erase old hexaloop entries */
+ memset(&Hexaloops, 0, 361);
+ memset(&Hexaloop37, 0, sizeof(int)*40);
+ memset(&HexaloopdH, 0, sizeof(int)*40);
+ do {
+ buf = get_line(fp);
+ if (buf==NULL) break;
+ r = sscanf(buf,"%8s %d %d", &Hexaloops[9*i], &Hexaloop37[i], &HexaloopdH[i]);
+ strcat(Hexaloops, " ");
+ free(buf);
+ i++;
+ } while((r==3)&&(i<40));
+ return;
+}
+
+/*------------------------------------------------------------*/
+PRIVATE void rd_Triloop37(void)
+{
+ int i, r;
+ char *buf;
+
+ i=0;
+ /* erase old hexaloop entries */
+ memset(&Triloops, 0, 241);
+ memset(&Triloop37, 0, sizeof(int)*40);
+ memset(&TriloopdH, 0, sizeof(int)*40);
+ do {
+ buf = get_line(fp);
+ if (buf==NULL) break;
+ r = sscanf(buf,"%5s %d %d", &Triloops[6*i], &Triloop37[i], &TriloopdH[i]);
+ strcat(Triloops, " ");
+ free(buf);
+ i++;
+ } while((r==3)&&(i<40));
+ return;
+}
+
+/*------------------------------------------------------------*/
+
+
+PRIVATE void ignore_comment(char * line)
+{
+ /* excise C style comments */
+ /* only one comment per line, no multiline comments */
+ char *cp1, *cp2;
+
+ if ((cp1=strstr(line, "/*"))) {
+ cp2 = strstr(cp1, "*/");
+ if (cp2==NULL)
+ nrerror("unclosed comment in parameter file");
+ /* can't use strcpy for overlapping strings */
+ for (cp2+=2; *cp2!='\0'; cp2++, cp1++)
+ *cp1 = *cp2;
+ *cp1 = '\0';
+ }
+
+ return;
+}
+/*------------------------------------------------------------*/
+
+PUBLIC char *settype(enum parset s){
+ switch(s){
+ case S: return "stack";
+ case S_H: return "stack_enthalpies";
+ case HP: return "hairpin";
+ case HP_H: return "hairpin_enthalpies";
+ case B: return "bulge";
+ case B_H: return "bulge_enthalpies";
+ case IL: return "interior";
+ case IL_H: return "interior_enthalpies";
+ case MME: return "mismatch_exterior";
+ case MME_H: return "mismatch_exterior_enthalpies";
+ case MMH: return "mismatch_hairpin";
+ case MMH_H: return "mismatch_hairpin_enthalpies";
+ case MMI: return "mismatch_interior";
+ case MMI_H: return "mismatch_interior_enthalpies";
+ case MMI1N: return "mismatch_interior_1n";
+ case MMI1N_H: return "mismatch_interior_1n_enthalpies";
+ case MMI23: return "mismatch_interior_23";
+ case MMI23_H: return "mismatch_interior_23_enthalpies";
+ case MMM: return "mismatch_multi";
+ case MMM_H: return "mismatch_multi_enthalpies";
+ case D5: return "dangle5";
+ case D5_H: return "dangle5_enthalpies";
+ case D3: return "dangle3";
+ case D3_H: return "dangle3_enthalpies";
+ case INT11: return "int11";
+ case INT11_H: return "int11_enthalpies";
+ case INT21: return "int21";
+ case INT21_H: return "int21_enthalpies";
+ case INT22: return "int22";
+ case INT22_H: return "int22_enthalpies";
+ case ML: return "ML_params";
+ case NIN: return "NINIO";
+ case TRI: return "Triloops";
+ case TL: return "Tetraloops";
+ case HEX: return "Hexaloops";
+ case QUIT: return "END";
+ case MISC: return "Misc";
+ default: nrerror("\nThe answer is: 42\n");
+ }
+ return "";
+}
+/*------------------------------------------------------------*/
+
+PUBLIC enum parset gettype(const char *ident){
+ if (strcmp(ident,"stack") == 0) return S;
+ else if (strcmp(ident,"stack_enthalpies") == 0) return S_H;
+ else if (strcmp(ident,"hairpin") == 0) return HP;
+ else if (strcmp(ident,"hairpin_enthalpies") == 0) return HP_H;
+ else if (strcmp(ident,"bulge") == 0) return B;
+ else if (strcmp(ident,"bulge_enthalpies") == 0) return B_H;
+ else if (strcmp(ident,"interior") == 0) return IL;
+ else if (strcmp(ident,"interior_enthalpies") == 0) return IL_H;
+ else if (strcmp(ident,"mismatch_exterior") == 0) return MME;
+ else if (strcmp(ident,"mismatch_exterior_enthalpies") == 0) return MME_H;
+ else if (strcmp(ident,"mismatch_hairpin") == 0) return MMH;
+ else if (strcmp(ident,"mismatch_hairpin_enthalpies") == 0) return MMH_H;
+ else if (strcmp(ident,"mismatch_interior") == 0) return MMI;
+ else if (strcmp(ident,"mismatch_interior_enthalpies") == 0) return MMI_H;
+ else if (strcmp(ident,"mismatch_interior_1n") == 0) return MMI1N;
+ else if (strcmp(ident,"mismatch_interior_1n_enthalpies") == 0) return MMI1N_H;
+ else if (strcmp(ident,"mismatch_interior_23") == 0) return MMI23;
+ else if (strcmp(ident,"mismatch_interior_23_enthalpies") == 0) return MMI23_H;
+ else if (strcmp(ident,"mismatch_multi") == 0) return MMM;
+ else if (strcmp(ident,"mismatch_multi_enthalpies") == 0) return MMM_H;
+ else if (strcmp(ident,"int11") == 0) return INT11;
+ else if (strcmp(ident,"int11_enthalpies") == 0) return INT11_H;
+ else if (strcmp(ident,"int21") == 0) return INT21;
+ else if (strcmp(ident,"int21_enthalpies") == 0) return INT21_H;
+ else if (strcmp(ident,"int22") == 0) return INT22;
+ else if (strcmp(ident,"int22_enthalpies") == 0) return INT22_H;
+ else if (strcmp(ident,"dangle5")== 0) return D5;
+ else if (strcmp(ident,"dangle5_enthalpies")== 0) return D5_H;
+ else if (strcmp(ident,"dangle3")== 0) return D3;
+ else if (strcmp(ident,"dangle3_enthalpies")== 0) return D3_H;
+ else if (strcmp(ident,"ML_params")== 0) return ML;
+ else if (strcmp(ident,"NINIO") == 0) return NIN;
+ else if (strcmp(ident,"Triloops") == 0) return TRI;
+ else if (strcmp(ident,"Tetraloops") == 0) return TL;
+ else if (strcmp(ident,"Hexaloops") == 0) return HEX;
+ else if (strcmp(ident,"Misc") == 0) return MISC;
+ else if (strcmp(ident,"END") == 0) return QUIT;
+ else return UNKNOWN;
+}
+
+/*---------------------------------------------------------------*/
+
+PUBLIC void write_parameter_file(const char fname[]){
+ FILE *outfp;
+ int c;
+ char *pnames[] = {"NP", "CG", "GC", "GU", "UG", "AU", "UA", " @"};
+ char bnames[] = "@ACGU";
+ outfp = fopen(fname, "w");
+ if (!outfp) {
+ fprintf(stderr, "can't open file %s\n", fname);
+ exit(1);
+ }
+ fprintf(outfp,"## RNAfold parameter file v2.0\n");
+
+ fprintf(outfp,"\n# %s\n", settype(S));
+ fprintf(outfp,"/* CG GC GU UG AU UA @ */\n");
+ for (c=1; c<NBPAIRS+1; c++)
+ display_array(stack37[c]+1,NBPAIRS,NBPAIRS, outfp);
+
+ fprintf(outfp,"\n# %s\n", settype(S_H));
+ fprintf(outfp,"/* CG GC GU UG AU UA @ */\n");
+ for (c=1; c<NBPAIRS+1; c++)
+ display_array(stackdH[c]+1,NBPAIRS,NBPAIRS, outfp);
+
+ fprintf(outfp,"\n# %s\n", settype(MMH));
+ { int i,k;
+ for (k=1; k<NBPAIRS+1; k++)
+ for (i=0; i<5; i++)
+ display_array(mismatchH37[k][i],5,5, outfp);
+ }
+
+ fprintf(outfp,"\n# %s\n", settype(MMH_H));
+ { int i,k;
+ for (k=1; k<NBPAIRS+1; k++)
+ for (i=0; i<5; i++)
+ display_array(mismatchHdH[k][i],5,5, outfp);
+
+ }
+
+ fprintf(outfp,"\n# %s\n", settype(MMI));
+ { int i,k;
+ for (k=1; k<NBPAIRS+1; k++)
+ for (i=0; i<5; i++)
+ display_array(mismatchI37[k][i],5,5, outfp);
+ }
+
+ fprintf(outfp,"\n# %s\n", settype(MMI_H));
+ { int i,k;
+ for (k=1; k<NBPAIRS+1; k++)
+ for (i=0; i<5; i++)
+ display_array(mismatchIdH[k][i],5,5, outfp);
+ }
+
+ fprintf(outfp,"\n# %s\n", settype(MMI1N));
+ { int i,k;
+ for (k=1; k<NBPAIRS+1; k++)
+ for (i=0; i<5; i++)
+ display_array(mismatch1nI37[k][i],5,5, outfp);
+ }
+
+ fprintf(outfp,"\n# %s\n", settype(MMI1N_H));
+ { int i,k;
+ for (k=1; k<NBPAIRS+1; k++)
+ for (i=0; i<5; i++)
+ display_array(mismatch1nIdH[k][i],5,5, outfp);
+ }
+
+ fprintf(outfp,"\n# %s\n", settype(MMI23));
+ { int i,k;
+ for (k=1; k<NBPAIRS+1; k++)
+ for (i=0; i<5; i++)
+ display_array(mismatch23I37[k][i],5,5, outfp);
+ }
+
+ fprintf(outfp,"\n# %s\n", settype(MMI23_H));
+ { int i,k;
+ for (k=1; k<NBPAIRS+1; k++)
+ for (i=0; i<5; i++)
+ display_array(mismatch23IdH[k][i],5,5, outfp);
+ }
+
+ fprintf(outfp,"\n# %s\n", settype(MMM));
+ { int i,k;
+ for (k=1; k<NBPAIRS+1; k++)
+ for (i=0; i<5; i++)
+ display_array(mismatchM37[k][i],5,5, outfp);
+ }
+
+ fprintf(outfp,"\n# %s\n", settype(MMM_H));
+ { int i,k;
+ for (k=1; k<NBPAIRS+1; k++)
+ for (i=0; i<5; i++)
+ display_array(mismatchMdH[k][i],5,5, outfp);
+ }
+
+ fprintf(outfp,"\n# %s\n", settype(MME));
+ { int i,k;
+ for (k=1; k<NBPAIRS+1; k++)
+ for (i=0; i<5; i++)
+ display_array(mismatchExt37[k][i],5,5, outfp);
+ }
+
+ fprintf(outfp,"\n# %s\n", settype(MME_H));
+ { int i,k;
+ for (k=1; k<NBPAIRS+1; k++)
+ for (i=0; i<5; i++)
+ display_array(mismatchExtdH[k][i],5,5, outfp);
+ }
+
+ fprintf(outfp,"\n# %s\n", settype(D5));
+ fprintf(outfp,"/* @ A C G U */\n");
+ for (c=1; c<NBPAIRS+1; c++)
+ display_array(dangle5_37[c], 5, 5, outfp);
+
+ fprintf(outfp,"\n# %s\n", settype(D5_H));
+ fprintf(outfp,"/* @ A C G U */\n");
+ for (c=1; c<NBPAIRS+1; c++)
+ display_array(dangle5_dH[c], 5, 5, outfp);
+
+ fprintf(outfp,"\n# %s\n", settype(D3));
+ fprintf(outfp,"/* @ A C G U */\n");
+ for (c=1; c<NBPAIRS+1; c++)
+ display_array(dangle3_37[c], 5, 5, outfp);
+
+ fprintf(outfp,"\n# %s\n", settype(D3_H));
+ fprintf(outfp,"/* @ A C G U */\n");
+ for (c=1; c<NBPAIRS+1; c++)
+ display_array(dangle3_dH[c], 5, 5, outfp);
+
+
+ /* dont print "no pair" entries for interior loop arrays */
+ fprintf(outfp,"\n# %s\n", settype(INT11));
+ { int i,k,l;
+ for (k=1; k<NBPAIRS+1; k++)
+ for (l=1; l<NBPAIRS+1; l++) {
+ fprintf(outfp, "/* %2s..%2s */\n", pnames[k], pnames[l]);
+ for (i=0; i<5; i++)
+ display_array(int11_37[k][l][i], 5, 5, outfp);
+ }
+ }
+
+ fprintf(outfp,"\n# %s\n", settype(INT11_H));
+ { int i,k,l;
+ for (k=1; k<NBPAIRS+1; k++)
+ for (l=1; l<NBPAIRS+1; l++) {
+ fprintf(outfp, "/* %2s..%2s */\n", pnames[k], pnames[l]);
+ for (i=0; i<5; i++)
+ display_array(int11_dH[k][l][i],5,5, outfp);
+ }
+ }
+
+ fprintf(outfp,"\n# %s\n", settype(INT21));
+ { int p1, p2, i, j;
+ for (p1=1; p1<NBPAIRS+1; p1++)
+ for (p2=1; p2<NBPAIRS+1; p2++)
+ for (i=0; i<5; i++) {
+ fprintf(outfp, "/* %2s.%c..%2s */\n",
+ pnames[p1], bnames[i], pnames[p2]);
+ for (j=0; j<5; j++)
+ display_array(int21_37[p1][p2][i][j],5,5, outfp);
+ }
+ }
+
+ fprintf(outfp,"\n# %s\n", settype(INT21_H));
+ { int p1, p2, i, j;
+ for (p1=1; p1<NBPAIRS+1; p1++)
+ for (p2=1; p2<NBPAIRS+1; p2++)
+ for (i=0; i<5; i++) {
+ fprintf(outfp, "/* %2s.%c..%2s */\n",
+ pnames[p1], bnames[i], pnames[p2]);
+ for (j=0; j<5; j++)
+ display_array(int21_dH[p1][p2][i][j],5,5, outfp);
+ }
+ }
+
+ fprintf(outfp,"\n# %s\n", settype(INT22));
+ { int p1, p2, i, j, k;
+ for (p1=1; p1<NBPAIRS; p1++)
+ for (p2=1; p2<NBPAIRS; p2++)
+ for (i=1; i<5; i++)
+ for (j=1; j<5; j++) {
+ fprintf(outfp, "/* %2s.%c%c..%2s */\n",
+ pnames[p1], bnames[i], bnames[j], pnames[p2]);
+ for (k=1; k<5; k++)
+ display_array(int22_37[p1][p2][i][j][k]+1,4,5, outfp);
+ }
+ }
+
+ fprintf(outfp,"\n# %s\n", settype(INT22_H));
+ { int p1, p2, i, j, k;
+ for (p1=1; p1<NBPAIRS; p1++)
+ for (p2=1; p2<NBPAIRS; p2++)
+ for (i=1; i<5; i++)
+ for (j=1; j<5; j++) {
+ fprintf(outfp, "/* %2s.%c%c..%2s */\n",
+ pnames[p1], bnames[i], bnames[j], pnames[p2]);
+ for (k=1; k<5; k++)
+ display_array(int22_dH[p1][p2][i][j][k]+1,4,5, outfp);
+ }
+ }
+
+ fprintf(outfp,"\n# %s\n", settype(HP));
+ display_array(hairpin37, 31, 10, outfp);
+
+ fprintf(outfp,"\n# %s\n", settype(HP_H));
+ display_array(hairpindH, 31, 10, outfp);
+
+ fprintf(outfp,"\n# %s\n", settype(B));
+ display_array(bulge37, 31, 10, outfp);
+
+ fprintf(outfp,"\n# %s\n", settype(B_H));
+ display_array(bulgedH, 31, 10, outfp);
+
+ fprintf(outfp,"\n# %s\n", settype(IL));
+ display_array(internal_loop37, 31, 10, outfp);
+
+ fprintf(outfp,"\n# %s\n", settype(IL_H));
+ display_array(internal_loopdH, 31, 10, outfp);
+
+ fprintf(outfp,"\n# %s\n", settype(ML));
+ fprintf(outfp,"/* F = cu*n_unpaired + cc + ci*loop_degree (+TermAU) */\n");
+ fprintf(outfp,"/*\t cu\t cu_dH\t cc\t cc_dH\t ci\t ci_dH */\n");
+ 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);
+
+ fprintf(outfp,"\n# %s\n", settype(NIN));
+ fprintf(outfp,"/* Ninio = MIN(max, m*|n1-n2| */\n"
+ "/*\t m\t m_dH max */\n"
+ "\t%6d\t%6d\t%6d\n", ninio37, niniodH, MAX_NINIO);
+
+ fprintf(outfp,"\n# %s\n", settype(MISC));
+ fprintf(outfp,"/* all parameters are pairs of 'energy enthalpy' */\n");
+ fprintf(outfp,"/* DuplexInit TerminalAU LXC */\n");
+ fprintf(outfp," %6d %6d %6d %6d %3.6f %6d\n", DuplexInit37, DuplexInitdH, TerminalAU37, TerminalAUdH, lxc37, 0);
+
+ fprintf(outfp,"\n# %s\n", settype(HEX));
+ for (c=0; c< strlen(Hexaloops)/9; c++)
+ fprintf(outfp,"\t%.8s %6d %6d\n", Hexaloops+c*9, Hexaloop37[c], HexaloopdH[c]);
+
+ fprintf(outfp,"\n# %s\n", settype(TL));
+ for (c=0; c< strlen(Tetraloops)/7; c++)
+ fprintf(outfp,"\t%.6s %6d %6d\n", Tetraloops+c*7, Tetraloop37[c], TetraloopdH[c]);
+
+ fprintf(outfp,"\n# %s\n", settype(TRI));
+ for (c=0; c< strlen(Triloops)/6; c++)
+ fprintf(outfp,"\t%.5s %6d %6d\n", Triloops+c*6, Triloop37[c], TriloopdH[c]);
+
+ fprintf(outfp,"\n# %s\n", settype(QUIT));
+ fclose(outfp);
+}
+
+PRIVATE void check_symmetry(void) {
+ int i,j,k,l;
+
+ for (i=0; i<=NBPAIRS; i++)
+ for (j=0; j<=NBPAIRS; j++)
+ if (stack37[i][j] != stack37[j][i])
+ fprintf(stderr, "WARNING: stacking energies not symmetric\n");
+
+ for (i=0; i<=NBPAIRS; i++)
+ for (j=0; j<=NBPAIRS; j++)
+ if (stackdH[i][j] != stackdH[j][i])
+ fprintf(stderr, "WARNING: stacking enthalpies not symmetric\n");
+
+
+ /* interior 1x1 loops */
+ for (i=0; i<=NBPAIRS; i++)
+ for (j=0; j<=NBPAIRS; j++)
+ for (k=0; k<5; k++)
+ for (l=0; l<5; l++)
+ if (int11_37[i][j][k][l] != int11_37[j][i][l][k])
+ 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]);
+
+ for (i=0; i<=NBPAIRS; i++)
+ for (j=0; j<=NBPAIRS; j++)
+ for (k=0; k<5; k++)
+ for (l=0; l<5; l++)
+ if (int11_dH[i][j][k][l] != int11_dH[j][i][l][k])
+ fprintf(stderr, "WARNING: int11 enthalpies not symmetric\n");
+
+ /* interior 2x2 loops */
+ for (i=0; i<=NBPAIRS; i++)
+ for (j=0; j<=NBPAIRS; j++)
+ for (k=0; k<5; k++)
+ for (l=0; l<5; l++) {
+ int m,n;
+ for (m=0; m<5; m++)
+ for (n=0; n<5; n++)
+ if (int22_37[i][j][k][l][m][n] != int22_37[j][i][m][n][k][l])
+ fprintf(stderr, "WARNING: int22 energies not symmetric\n");
+ }
+
+ for (i=0; i<=NBPAIRS; i++)
+ for (j=0; j<=NBPAIRS; j++)
+ for (k=0; k<5; k++)
+ for (l=0; l<5; l++) {
+ int m,n;
+ for (m=0; m<5; m++)
+ for (n=0; n<5; n++)
+ if (int22_dH[i][j][k][l][m][n] != int22_dH[j][i][m][n][k][l])
+ fprintf(stderr, "WARNING: int22 enthalpies not symmetric: %d %d %d %d %d %d\n", i,j,k,l,m,n);
+ }
+}
+
+/* update nonstandard nucleotide/basepair involved contributions for int22 */
+PRIVATE void update_nst(int array[NBPAIRS+1][NBPAIRS+1][5][5][5][5]){
+ int i, j, k, l, m, n;
+ int max, max2, max3, max4, max5, max6;
+
+ /* get maxima for one nonstandard nucleotide */
+ for (i=1; i<NBPAIRS; i++){
+ for (j=1; j<NBPAIRS; j++){
+ for (k=1; k<5; k++){
+ for (l=1; l<5; l++){
+ for (m=1; m<5; m++){
+ max = max2 = max3 = max4 = -INF; /* max of {CGAU} */
+ for(n=1;n<5;n++){
+ max = MAX2(max, array[i][j][k][l][m][n]);
+ max2 = MAX2(max2, array[i][j][k][l][n][m]);
+ max3 = MAX2(max3, array[i][j][k][n][l][m]);
+ max4 = MAX2(max4, array[i][j][n][k][l][m]);
+ }
+ array[i][j][k][l][m][0] = max;
+ array[i][j][k][l][0][m] = max2;
+ array[i][j][k][0][l][m] = max3;
+ array[i][j][0][k][l][m] = max4;
+ }
+ }
+ }
+ }
+ }
+ /* get maxima for two nonstandard nucleotides */
+ for (i=1; i<NBPAIRS; i++){
+ for (j=1; j<NBPAIRS; j++){
+ for (k=1; k<5; k++){
+ for (l=1; l<5; l++){
+ max = max2 = max3 = max4 = max5 = max6 = -INF; /* max of {CGAU} */
+ for (m=1; m<5; m++){
+ max = MAX2(max, array[i][j][k][l][m][0]);
+ max2 = MAX2(max2, array[i][j][k][m][0][l]);
+ max3 = MAX2(max3, array[i][j][m][0][k][l]);
+ max4 = MAX2(max4, array[i][j][0][k][l][m]);
+ max5 = MAX2(max5, array[i][j][0][k][m][l]);
+ max6 = MAX2(max6, array[i][j][k][0][l][m]);
+ }
+ array[i][j][k][l][0][0] = max;
+ array[i][j][k][0][0][l] = max2;
+ array[i][j][0][0][k][l] = max3;
+ array[i][j][k][0][l][0] = max6;
+ array[i][j][0][k][0][l] = max5;
+ array[i][j][0][k][l][0] = max4;
+ }
+ }
+ }
+ }
+ /* get maxima for three nonstandard nucleotides */
+ for (i=1; i<NBPAIRS; i++){
+ for (j=1; j<NBPAIRS; j++){
+ for (k=1; k<5; k++){
+ max = max2 = max3 = max4 = -INF; /* max of {CGAU} */
+ for (l=1; l<5; l++){
+ /* should be arbitrary where index l resides in last 3 possible locations */
+ max = MAX2(max, array[i][j][k][l][0][0]);
+ max2 = MAX2(max2, array[i][j][0][k][l][0]);
+ max3 = MAX2(max3, array[i][j][0][0][k][l]);
+ max4 = MAX2(max4, array[i][j][0][0][l][k]);
+ }
+ array[i][j][k][0][0][0] = max;
+ array[i][j][0][k][0][0] = max2;
+ array[i][j][0][0][k][0] = max3;
+ array[i][j][0][0][0][k] = max4;
+ }
+ }
+ }
+ /* get maxima for 4 nonstandard nucleotides */
+ for (i=1; i<NBPAIRS; i++){
+ for (j=1; j<NBPAIRS; j++){
+ max = -INF; /* max of {CGAU} */
+ for (k=1; k<5; k++){
+ max = MAX2(max, array[i][j][k][0][0][0]);
+ }
+ array[i][j][0][0][0][0] = max;
+ }
+ }
+
+ /* now compute contributions for nonstandard base pairs ... */
+ /* first, 1 nonstandard bp */
+ for (i=1; i<NBPAIRS; i++){
+ for (k=0; k<5; k++){
+ for (l=0; l<5; l++){
+ for (m=0; m<5; m++){
+ for(n=0;n<5;n++){
+ max = max2 = -INF;
+ for(j=1;j<NBPAIRS;j++){
+ max = MAX2(max, array[i][j][k][l][m][n]);
+ max2 = MAX2(max2, array[j][i][k][l][m][n]);
+ }
+ array[i][NBPAIRS][k][l][m][n] = max;
+ array[NBPAIRS][i][k][l][m][n] = max2;
+ }
+ }
+ }
+ }
+ }
+
+ /* now 2 nst base pairs */
+ for (k=0; k<5; k++){
+ for (l=0; l<5; l++){
+ for (m=0; m<5; m++){
+ for(n=0;n<5;n++){
+ max = -INF;
+ for(j=1;j<NBPAIRS;j++){
+ max = MAX2(max, array[NBPAIRS][j][k][l][m][n]);
+ }
+ array[NBPAIRS][NBPAIRS][k][l][m][n] = max;
+ }
+ }
+ }
+ }
+
+}
git://source.jalview.org / jabaws.git / blobdiff