/*
                  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"

static char rcsid[] = "$Id: read_epars.c,v 1.10 2004/12/10 16:32:35 ivo Exp $";

#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;
        }
      }
    }
  }

}