binaries/src/ViennaRNA/lib/read_epars.c

   1 /*
   2                   read energy parameters from a file
   3
   4                       Stephan Kopp, Ivo Hofacker
   5                           Vienna RNA Package
   6 */
   7
   8
   9 #include <stdio.h>
  10 #include <stdlib.h>
  11 #include <ctype.h>
  12 #include <string.h>
  13 #include <math.h>
  14 #include <stdarg.h>
  15 #include "utils.h"
  16 #include "energy_const.h"
  17 #include "energy_par.h"
  18 #include "read_epars.h"
  19
  20 #define PUBLIC
  21 #define PRIVATE   static
  22 #define PARSET 20
  23
  24 #define DEF -50
  25 #define NST 0
  26
  27 PRIVATE FILE *fp;
  28
  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]);
  34
  35 /**
  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
  40 **/
  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,
  44                       int dim1, int dim2,
  45                       int shift1, int shift2);
  46 PRIVATE void  rd_2dim_slice(int *array,
  47                       int dim1, int dim2,
  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);
  81
  82 /*------------------------------------------------------------*/
  83 PUBLIC void read_parameter_file(const char fname[]){
  84   char        *line, ident[256];
  85   enum parset type;
  86   int         r;
  87
  88   if (!(fp=fopen(fname,"r"))) {
  89     fprintf(stderr,
  90             "\nread_parameter_file:\n"
  91             "\t\tcan't open file %s\n"
  92             "\t\tusing default parameters instead.\n", fname);
  93     return;
  94   }
  95
  96   if (!(line = get_line(fp))) {
  97     fprintf(stderr," File %s is inproper.\n", fname);
  98     fclose(fp);
  99     return;
 100   }
 101
 102   if (strncmp(line,"## RNAfold parameter file v2.0",30)!=0) {
 103     fprintf(stderr,
 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");
 107   }
 108   free(line);
 109
 110   while((line=get_line(fp))) {
 111
 112     r = sscanf(line, "# %255s", ident);
 113     if (r==1) {
 114       type = gettype(ident);
 115       switch (type){
 116         case QUIT:    break;
 117         case S:       rd_2dim(&(stack37[0][0]), NBPAIRS+1, NBPAIRS+1, 1, 1);
 118                       break;
 119         case S_H:     rd_2dim(&(stackdH[0][0]), NBPAIRS+1, NBPAIRS+1, 1, 1);
 120                       break;
 121         case HP:      rd_1dim(&(hairpin37[0]), 31, 0);
 122                       break;
 123         case HP_H:    rd_1dim(&(hairpindH[0]), 31, 0);
 124                       break;
 125         case B:       rd_1dim(&(bulge37[0]), 31, 0);
 126                       break;
 127         case B_H:     rd_1dim(&(bulgedH[0]), 31, 0);
 128                       break;
 129         case IL:      rd_1dim(&(internal_loop37[0]), 31, 0);
 130                       break;
 131         case IL_H:    rd_1dim(&(internal_loopdH[0]), 31, 0);
 132                       break;
 133         case MME:     rd_3dim(&(mismatchExt37[0][0][0]),
 134                           NBPAIRS+1, 5, 5,
 135                           1, 0, 0);
 136                       break;
 137         case MME_H:   rd_3dim(&(mismatchExtdH[0][0][0]),
 138                           NBPAIRS+1, 5, 5,
 139                           1, 0, 0);
 140                       break;
 141         case MMH:     rd_3dim(&(mismatchH37[0][0][0]),
 142                           NBPAIRS+1, 5, 5,
 143                           1, 0, 0);
 144                       break;
 145         case MMH_H:   rd_3dim(&(mismatchHdH[0][0][0]),
 146                           NBPAIRS+1, 5, 5,
 147                           1, 0, 0);
 148                       break;
 149         case MMI:     rd_3dim(&(mismatchI37[0][0][0]),
 150                           NBPAIRS+1, 5, 5,
 151                           1, 0, 0);
 152                       break;
 153         case MMI_H:   rd_3dim(&(mismatchIdH[0][0][0]),
 154                           NBPAIRS+1, 5, 5,
 155                           1, 0, 0);
 156                       break;
 157         case MMI1N:   rd_3dim(&(mismatch1nI37[0][0][0]),
 158                           NBPAIRS+1, 5, 5,
 159                           1, 0, 0);
 160                       break;
 161         case MMI1N_H: rd_3dim(&(mismatch1nIdH[0][0][0]),
 162                           NBPAIRS+1, 5, 5,
 163                           1, 0, 0);
 164                       break;
 165         case MMI23:   rd_3dim(&(mismatch23I37[0][0][0]),
 166                           NBPAIRS+1, 5, 5,
 167                           1, 0, 0);
 168                       break;
 169         case MMI23_H: rd_3dim(&(mismatch23IdH[0][0][0]),
 170                           NBPAIRS+1, 5, 5,
 171                           1, 0, 0);
 172                       break;
 173         case MMM:     rd_3dim(&(mismatchM37[0][0][0]),
 174                           NBPAIRS+1, 5, 5,
 175                           1, 0, 0);
 176                       break;
 177         case MMM_H:   rd_3dim(&(mismatchMdH[0][0][0]),
 178                           NBPAIRS+1, 5, 5,
 179                           1, 0, 0);
 180                       break;
 181         case INT11:   rd_4dim(&(int11_37[0][0][0][0]),
 182                           NBPAIRS+1, NBPAIRS+1, 5, 5,
 183                           1, 1, 0, 0);
 184                       break;
 185         case INT11_H: rd_4dim(&(int11_dH[0][0][0][0]),
 186                           NBPAIRS+1, NBPAIRS+1, 5, 5,
 187                           1, 1, 0, 0);
 188                       break;
 189         case INT21:   rd_5dim(&(int21_37[0][0][0][0][0]),
 190                           NBPAIRS+1, NBPAIRS+1, 5, 5, 5,
 191                           1, 1, 0, 0, 0);
 192                       break;
 193         case INT21_H: rd_5dim(&(int21_dH[0][0][0][0][0]),
 194                           NBPAIRS+1, NBPAIRS+1, 5, 5, 5,
 195                           1, 1, 0, 0, 0);
 196                       break;
 197         case INT22:   rd_6dim_slice(&(int22_37[0][0][0][0][0][0]),
 198                           NBPAIRS+1, NBPAIRS+1, 5, 5, 5, 5,
 199                           1, 1, 1, 1, 1, 1,
 200                           1, 1, 0, 0, 0, 0);
 201                       update_nst(int22_37);
 202                       break;
 203         case INT22_H: rd_6dim_slice(&(int22_dH[0][0][0][0][0][0]),
 204                           NBPAIRS+1, NBPAIRS+1, 5, 5, 5, 5,
 205                           1, 1, 1, 1, 1, 1,
 206                           1, 1, 0, 0, 0, 0);
 207                       update_nst(int22_dH);
 208                       break;
 209         case D5:      rd_2dim(&(dangle5_37[0][0]), NBPAIRS+1, 5, 1, 0);
 210                       break;
 211         case D5_H:    rd_2dim(&(dangle5_dH[0][0]), NBPAIRS+1, 5, 1, 0);
 212                       break;
 213         case D3:      rd_2dim(&(dangle3_37[0][0]), NBPAIRS+1, 5, 1, 0);
 214                       break;
 215         case D3_H:    rd_2dim(&(dangle3_dH[0][0]), NBPAIRS+1, 5, 1, 0);
 216                       break;
 217         case ML:      {
 218                         int values[6];
 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];
 226                       }
 227                       break;
 228         case NIN:     {
 229                         int values[3];
 230                         rd_1dim(&values[0], 3, 0);
 231                         ninio37 = values[0];
 232                         niniodH = values[1];
 233                         MAX_NINIO  = values[2];
 234                       }
 235                       break;
 236         case MISC:    {
 237                         int values[4];
 238                         rd_1dim(&values[0], 4, 0);
 239                         DuplexInit37 = values[0];
 240                         DuplexInitdH = values[1];
 241                         TerminalAU37 = values[2];
 242                         TerminalAUdH = values[3];
 243                       }
 244                       break;
 245         case TL:      rd_Tetraloop37();
 246                       break;
 247         case TRI:     rd_Triloop37();
 248                       break;
 249         case HEX:     rd_Hexaloop37();
 250                       break;
 251         default:      /* do nothing but complain */
 252                       fprintf(stderr,"read_epars: Unknown field identifier in `%s'\n", line);
 253       }
 254     } /* else ignore line */
 255     free(line);
 256   }
 257   fclose(fp);
 258
 259   check_symmetry();
 260   return;
 261 }
 262
 263 /*------------------------------------------------------------*/
 264
 265 PRIVATE void display_array(int *p, int size, int nl, FILE *fp)
 266 {
 267   int i;
 268
 269   for (i=1; i<=size; i++, p++) {
 270     switch(*p)
 271       {
 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;
 276       }
 277     if ((i%nl)==0) fprintf(fp,"\n");
 278   }
 279   if (size%nl) fprintf(fp,"\n");
 280   return;
 281 }
 282
 283 /*------------------------------------------------------------*/
 284
 285 PRIVATE char *get_array1(int *arr, int size)
 286 {
 287   int    i, p, pos, pp, r, last;
 288   char  *line, buf[16];
 289
 290
 291   i = last = 0;
 292   while( i<size ) {
 293     line = get_line(fp);
 294     if (!line) nrerror("unexpected end of file in get_array1");
 295     ignore_comment(line);
 296     pos=0;
 297     while ((i<size)&&(sscanf(line+pos,"%15s%n", buf, &pp)==1)) {
 298       pos += pp;
 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)));
 303       }
 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;
 307       else {
 308         r=sscanf(buf,"%d", &p);
 309         if (r!=1) {
 310           return line+pos;
 311           fprintf(stderr, "can't interpret `%s' in get_array1\n", buf);
 312           exit(1);
 313         }
 314         last = i;
 315       }
 316       arr[i++]=p;
 317     }
 318     free(line);
 319   }
 320
 321   return NULL;
 322 }
 323
 324 PRIVATE void rd_1dim(int *array, int dim, int shift){
 325   rd_1dim_slice(array, dim, shift, 0);
 326 }
 327
 328 PRIVATE void rd_1dim_slice(int *array, int dim, int shift, int post){
 329   char *cp;
 330   cp   = get_array1(array+shift, dim-shift-post);
 331
 332   if (cp) {
 333     fprintf(stderr,"\nrd_1dim: %s\n", cp);
 334     exit(1);
 335   }
 336   return;
 337 }
 338
 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);
 341 }
 342
 343 PRIVATE void  rd_2dim_slice(int *array,
 344                       int dim1, int dim2,
 345                       int shift1, int shift2,
 346                       int post1, int post2){
 347   int i;
 348   int delta_pre   = shift1 + shift2;
 349   int delta_post  = post1 + post2;
 350
 351   if(delta_pre + delta_post == 0){
 352     rd_1dim(array, dim1 * dim2, 0);
 353     return;
 354   }
 355   for (i=shift1; i<dim1 - post1; i++)
 356     rd_1dim_slice(array + (i*dim2), dim2, shift2, post2);
 357   return;
 358 }
 359
 360 PRIVATE void  rd_3dim(int *array, int dim1, int dim2, int dim3, int shift1, int shift2, int shift3){
 361   rd_3dim_slice(array,
 362                 dim1, dim2, dim3,
 363                 shift1, shift2, shift3,
 364                 0, 0, 0);
 365 }
 366
 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){
 371   int    i;
 372   int delta_pre   = shift1 + shift2 + shift3;
 373   int delta_post  = post1 + post2 + post3;
 374
 375   if(delta_pre + delta_post == 0){
 376     rd_1dim(array, dim1 * dim2 * dim3, 0);
 377     return;
 378   }
 379   for (i=shift1; i<dim1 - post1; i++){
 380     rd_2dim_slice(array + (i * dim2 * dim3),
 381             dim2, dim3,
 382             shift2, shift3,
 383             post2, post3);
 384   }
 385   return;
 386 }
 387
 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){
 391   rd_4dim_slice(array,
 392                 dim1, dim2, dim3, dim4,
 393                 shift1, shift2, shift3, shift4,
 394                 0, 0, 0, 0);
 395 }
 396
 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){
 401   int i;
 402   int delta_pre   = shift1 + shift2 + shift3 + shift4;
 403   int delta_post  = post1 + post2 + post3 + post4;
 404
 405   if(delta_pre + delta_post == 0){
 406     rd_1dim(array, dim1 * dim2 * dim3 * dim4, 0);
 407     return;
 408   }
 409   for(i=shift1; i<dim1 - post1; i++){
 410     rd_3dim_slice(array + (i * dim2 * dim3 * dim4),
 411             dim2, dim3, dim4,
 412             shift2, shift3, shift4,
 413             post2, post3, post4);
 414   }
 415   return;
 416 }
 417
 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){
 421   rd_5dim_slice(array,
 422                 dim1, dim2, dim3, dim4, dim5,
 423                 shift1, shift2, shift3, shift4, shift5,
 424                 0, 0, 0, 0, 0);
 425 }
 426
 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){
 431   int i;
 432   int delta_pre   = shift1 + shift2 + shift3 + shift4 + shift5;
 433   int delta_post  = post1 + post2 + post3 + post4 + post5;
 434
 435   if(delta_pre + delta_post == 0){
 436     rd_1dim(array, dim1 * dim2 * dim3 * dim4 * dim5, 0);
 437     return;
 438   }
 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);
 444   return;
 445 }
 446
 447 /**
 448 *** \param dim1   The size of the first dimension
 449 *** \param shift1 The pre shift for the first dimension
 450 **/
 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){
 454   rd_6dim_slice(array,
 455                 dim1, dim2, dim3, dim4, dim5, dim6,
 456                 shift1, shift2, shift3, shift4, shift5, shift6,
 457                 0, 0, 0, 0, 0, 0);
 458 }
 459
 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){
 464   int i;
 465   int delta_pre   = shift1 + shift2 + shift3 + shift4 + shift5 + shift6;
 466   int delta_post  = post1 + post2 + post3 + post4 + post5 + post6;
 467
 468   if(delta_pre + delta_post == 0){
 469     rd_1dim(array, dim1 * dim2 * dim3 * dim4 * dim5 * dim6, 0);
 470     return;
 471   }
 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);
 477   return;
 478 }
 479
 480
 481 /*------------------------------------------------------------*/
 482 PRIVATE void  rd_Tetraloop37(void)
 483 {
 484   int    i, r;
 485   char   *buf;
 486
 487   i=0;
 488   /* erase old tetraloop entries */
 489   memset(&Tetraloops, 0, 281);
 490   memset(&Tetraloop37, 0, sizeof(int)*40);
 491   memset(&TetraloopdH, 0, sizeof(int)*40);
 492   do {
 493     buf = get_line(fp);
 494     if (buf==NULL) break;
 495     r = sscanf(buf,"%6s %d %d", &Tetraloops[7*i], &Tetraloop37[i], &TetraloopdH[i]);
 496     strcat(Tetraloops, " ");
 497     free(buf);
 498     i++;
 499   } while((r==3)&&(i<40));
 500   return;
 501 }
 502
 503 /*------------------------------------------------------------*/
 504 PRIVATE void  rd_Hexaloop37(void)
 505 {
 506   int    i, r;
 507   char   *buf;
 508
 509   i=0;
 510   /* erase old hexaloop entries */
 511   memset(&Hexaloops, 0, 361);
 512   memset(&Hexaloop37, 0, sizeof(int)*40);
 513   memset(&HexaloopdH, 0, sizeof(int)*40);
 514   do {
 515     buf = get_line(fp);
 516     if (buf==NULL) break;
 517     r = sscanf(buf,"%8s %d %d", &Hexaloops[9*i], &Hexaloop37[i], &HexaloopdH[i]);
 518     strcat(Hexaloops, " ");
 519     free(buf);
 520     i++;
 521   } while((r==3)&&(i<40));
 522   return;
 523 }
 524
 525 /*------------------------------------------------------------*/
 526 PRIVATE void  rd_Triloop37(void)
 527 {
 528   int    i, r;
 529   char   *buf;
 530
 531   i=0;
 532   /* erase old hexaloop entries */
 533   memset(&Triloops,   0, 241);
 534   memset(&Triloop37,  0, sizeof(int)*40);
 535   memset(&TriloopdH,  0, sizeof(int)*40);
 536   do {
 537     buf = get_line(fp);
 538     if (buf==NULL) break;
 539     r = sscanf(buf,"%5s %d %d", &Triloops[6*i], &Triloop37[i], &TriloopdH[i]);
 540     strcat(Triloops, " ");
 541     free(buf);
 542     i++;
 543   } while((r==3)&&(i<40));
 544   return;
 545 }
 546
 547 /*------------------------------------------------------------*/
 548
 549
 550 PRIVATE void ignore_comment(char * line)
 551 {
 552   /* excise C style comments */
 553   /* only one comment per line, no multiline comments */
 554   char *cp1, *cp2;
 555
 556   if ((cp1=strstr(line, "/*"))) {
 557     cp2 = strstr(cp1, "*/");
 558     if (cp2==NULL)
 559       nrerror("unclosed comment in parameter file");
 560     /* can't use strcpy for overlapping strings */
 561     for (cp2+=2; *cp2!='\0'; cp2++, cp1++)
 562       *cp1 = *cp2;
 563     *cp1 = '\0';
 564   }
 565
 566   return;
 567 }
 568 /*------------------------------------------------------------*/
 569
 570 PUBLIC char *settype(enum parset s){
 571   switch(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");
 610   }
 611   return "";
 612 }
 613 /*------------------------------------------------------------*/
 614
 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;
 653   else return UNKNOWN;
 654 }
 655
 656 /*---------------------------------------------------------------*/
 657
 658 PUBLIC void write_parameter_file(const char fname[]){
 659   FILE *outfp;
 660   int c;
 661   char *pnames[] = {"NP", "CG", "GC", "GU", "UG", "AU", "UA", " @"};
 662   char bnames[] = "@ACGU";
 663   outfp = fopen(fname, "w");
 664   if (!outfp) {
 665     fprintf(stderr, "can't open file %s\n", fname);
 666     exit(1);
 667   }
 668   fprintf(outfp,"## RNAfold parameter file v2.0\n");
 669
 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);
 674
 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);
 679
 680   fprintf(outfp,"\n# %s\n", settype(MMH));
 681   { int i,k;
 682   for (k=1; k<NBPAIRS+1; k++)
 683     for (i=0; i<5; i++)
 684       display_array(mismatchH37[k][i],5,5, outfp);
 685   }
 686
 687   fprintf(outfp,"\n# %s\n", settype(MMH_H));
 688   { int i,k;
 689   for (k=1; k<NBPAIRS+1; k++)
 690     for (i=0; i<5; i++)
 691       display_array(mismatchHdH[k][i],5,5, outfp);
 692
 693   }
 694
 695   fprintf(outfp,"\n# %s\n", settype(MMI));
 696   { int i,k;
 697   for (k=1; k<NBPAIRS+1; k++)
 698     for (i=0; i<5; i++)
 699       display_array(mismatchI37[k][i],5,5, outfp);
 700   }
 701
 702   fprintf(outfp,"\n# %s\n", settype(MMI_H));
 703   { int i,k;
 704   for (k=1; k<NBPAIRS+1; k++)
 705     for (i=0; i<5; i++)
 706       display_array(mismatchIdH[k][i],5,5, outfp);
 707   }
 708
 709   fprintf(outfp,"\n# %s\n", settype(MMI1N));
 710   { int i,k;
 711   for (k=1; k<NBPAIRS+1; k++)
 712     for (i=0; i<5; i++)
 713       display_array(mismatch1nI37[k][i],5,5, outfp);
 714   }
 715
 716   fprintf(outfp,"\n# %s\n", settype(MMI1N_H));
 717   { int i,k;
 718   for (k=1; k<NBPAIRS+1; k++)
 719     for (i=0; i<5; i++)
 720       display_array(mismatch1nIdH[k][i],5,5, outfp);
 721   }
 722
 723   fprintf(outfp,"\n# %s\n", settype(MMI23));
 724   { int i,k;
 725   for (k=1; k<NBPAIRS+1; k++)
 726     for (i=0; i<5; i++)
 727       display_array(mismatch23I37[k][i],5,5, outfp);
 728   }
 729
 730   fprintf(outfp,"\n# %s\n", settype(MMI23_H));
 731   { int i,k;
 732   for (k=1; k<NBPAIRS+1; k++)
 733     for (i=0; i<5; i++)
 734       display_array(mismatch23IdH[k][i],5,5, outfp);
 735   }
 736
 737   fprintf(outfp,"\n# %s\n", settype(MMM));
 738   { int i,k;
 739   for (k=1; k<NBPAIRS+1; k++)
 740     for (i=0; i<5; i++)
 741       display_array(mismatchM37[k][i],5,5, outfp);
 742   }
 743
 744   fprintf(outfp,"\n# %s\n", settype(MMM_H));
 745   { int i,k;
 746   for (k=1; k<NBPAIRS+1; k++)
 747     for (i=0; i<5; i++)
 748       display_array(mismatchMdH[k][i],5,5, outfp);
 749   }
 750
 751   fprintf(outfp,"\n# %s\n", settype(MME));
 752   { int i,k;
 753   for (k=1; k<NBPAIRS+1; k++)
 754     for (i=0; i<5; i++)
 755       display_array(mismatchExt37[k][i],5,5, outfp);
 756   }
 757
 758   fprintf(outfp,"\n# %s\n", settype(MME_H));
 759   { int i,k;
 760   for (k=1; k<NBPAIRS+1; k++)
 761     for (i=0; i<5; i++)
 762       display_array(mismatchExtdH[k][i],5,5, outfp);
 763   }
 764
 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);
 769
 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);
 774
 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);
 779
 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);
 784
 785
 786   /* dont print "no pair" entries for interior loop arrays */
 787   fprintf(outfp,"\n# %s\n", settype(INT11));
 788   { int i,k,l;
 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]);
 792       for (i=0; i<5; i++)
 793         display_array(int11_37[k][l][i], 5, 5, outfp);
 794     }
 795   }
 796
 797   fprintf(outfp,"\n# %s\n", settype(INT11_H));
 798   { int i,k,l;
 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]);
 802       for (i=0; i<5; i++)
 803         display_array(int11_dH[k][l][i],5,5, outfp);
 804     }
 805   }
 806
 807   fprintf(outfp,"\n# %s\n", settype(INT21));
 808   { int p1, p2, i, j;
 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]);
 814         for (j=0; j<5; j++)
 815           display_array(int21_37[p1][p2][i][j],5,5, outfp);
 816       }
 817   }
 818
 819   fprintf(outfp,"\n# %s\n", settype(INT21_H));
 820   { int p1, p2, i, j;
 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]);
 826         for (j=0; j<5; j++)
 827           display_array(int21_dH[p1][p2][i][j],5,5, outfp);
 828       }
 829   }
 830
 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++)
 835       for (i=1; i<5; i++)
 836         for (j=1; j<5; j++) {
 837           fprintf(outfp, "/* %2s.%c%c..%2s */\n",
 838                   pnames[p1], bnames[i], bnames[j], pnames[p2]);
 839           for (k=1; k<5; k++)
 840             display_array(int22_37[p1][p2][i][j][k]+1,4,5, outfp);
 841         }
 842   }
 843
 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++)
 848       for (i=1; i<5; i++)
 849         for (j=1; j<5; j++) {
 850           fprintf(outfp, "/* %2s.%c%c..%2s */\n",
 851                   pnames[p1], bnames[i], bnames[j], pnames[p2]);
 852           for (k=1; k<5; k++)
 853             display_array(int22_dH[p1][p2][i][j][k]+1,4,5, outfp);
 854         }
 855   }
 856
 857   fprintf(outfp,"\n# %s\n", settype(HP));
 858   display_array(hairpin37, 31, 10, outfp);
 859
 860   fprintf(outfp,"\n# %s\n", settype(HP_H));
 861   display_array(hairpindH, 31, 10, outfp);
 862
 863   fprintf(outfp,"\n# %s\n", settype(B));
 864   display_array(bulge37, 31, 10, outfp);
 865
 866   fprintf(outfp,"\n# %s\n", settype(B_H));
 867   display_array(bulgedH, 31, 10, outfp);
 868
 869   fprintf(outfp,"\n# %s\n", settype(IL));
 870   display_array(internal_loop37, 31, 10, outfp);
 871
 872   fprintf(outfp,"\n# %s\n", settype(IL_H));
 873   display_array(internal_loopdH, 31, 10, outfp);
 874
 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);
 879
 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);
 884
 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);
 889
 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]);
 893
 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]);
 897
 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]);
 901
 902   fprintf(outfp,"\n# %s\n", settype(QUIT));
 903   fclose(outfp);
 904 }
 905
 906 PRIVATE void check_symmetry(void) {
 907   int i,j,k,l;
 908
 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");
 913
 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");
 918
 919
 920   /* interior 1x1 loops */
 921   for (i=0; i<=NBPAIRS; i++)
 922     for (j=0; j<=NBPAIRS; j++)
 923       for (k=0; k<5; k++)
 924         for (l=0; l<5; l++)
 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]);
 927
 928   for (i=0; i<=NBPAIRS; i++)
 929     for (j=0; j<=NBPAIRS; j++)
 930       for (k=0; k<5; k++)
 931         for (l=0; l<5; l++)
 932           if (int11_dH[i][j][k][l] != int11_dH[j][i][l][k])
 933             fprintf(stderr, "WARNING: int11 enthalpies not symmetric\n");
 934
 935   /* interior 2x2 loops */
 936   for (i=0; i<=NBPAIRS; i++)
 937     for (j=0; j<=NBPAIRS; j++)
 938       for (k=0; k<5; k++)
 939         for (l=0; l<5; l++) {
 940           int m,n;
 941           for (m=0; m<5; m++)
 942             for (n=0; n<5; n++)
 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");
 945         }
 946
 947   for (i=0; i<=NBPAIRS; i++)
 948     for (j=0; j<=NBPAIRS; j++)
 949       for (k=0; k<5; k++)
 950         for (l=0; l<5; l++) {
 951           int m,n;
 952           for (m=0; m<5; m++)
 953             for (n=0; n<5; n++)
 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);
 956         }
 957 }
 958
 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;
 963
 964   /* get maxima for one nonstandard nucleotide */
 965   for (i=1; i<NBPAIRS; i++){
 966     for (j=1; j<NBPAIRS; j++){
 967       for (k=1; k<5; k++){
 968         for (l=1; l<5; l++){
 969           for (m=1; m<5; m++){
 970             max = max2 = max3 = max4 = -INF; /* max of {CGAU} */
 971             for(n=1;n<5;n++){
 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]);
 976             }
 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;
 981           }
 982         }
 983       }
 984     }
 985   }
 986   /* get maxima for two nonstandard nucleotides */
 987   for (i=1; i<NBPAIRS; i++){
 988     for (j=1; j<NBPAIRS; j++){
 989       for (k=1; k<5; k++){
 990         for (l=1; l<5; l++){
 991           max = max2 = max3 = max4 = max5 = max6 = -INF; /* max of {CGAU} */
 992           for (m=1; m<5; m++){
 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]);
 999           }
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;
1006         }
1007       }
1008     }
1009   }
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]);
1021         }
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;
1026       }
1027     }
1028   }
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]);
1035       }
1036       array[i][j][0][0][0][0] = max;
1037     }
1038   }
1039
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++){
1046           for(n=0;n<5;n++){
1047             max = max2 = -INF;
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]);
1051             }
1052             array[i][NBPAIRS][k][l][m][n] = max;
1053             array[NBPAIRS][i][k][l][m][n] = max2;
1054           }
1055         }
1056       }
1057     }
1058   }
1059
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++){
1064         for(n=0;n<5;n++){
1065           max = -INF;
1066           for(j=1;j<NBPAIRS;j++){
1067             max   = MAX2(max, array[NBPAIRS][j][k][l][m][n]);
1068           }
1069           array[NBPAIRS][NBPAIRS][k][l][m][n] = max;
1070         }
1071       }
1072     }
1073   }
1074
1075 }