pair_mat.h

#ifndef VIENNA_RNA_PACKAGE_PAIR_MAT_H
#define VIENNA_RNA_PACKAGE_PAIR_MAT_H

#include <ctype.h>
#include <ViennaRNA/utils/basic.h>
#include <ViennaRNA/fold_vars.h>

#define NBASES 8
/*@notnull@*/

#ifndef INLINE
# ifdef __GNUC__
#  define INLINE inline
# else
#  define INLINE
# endif
#endif

static const char Law_and_Order[]         = "_ACGUTXKI";
static int        BP_pair[NBASES][NBASES] =
  /* _  A  C  G  U  X  K  I */
{ { 0, 0, 0, 0, 0, 0, 0, 0 },
  { 0, 0, 0, 0, 5, 0, 0, 5 },
  { 0, 0, 0, 1, 0, 0, 0, 0 },
  { 0, 0, 2, 0, 3, 0, 0, 0 },
  { 0, 6, 0, 4, 0, 0, 0, 6 },
  { 0, 0, 0, 0, 0, 0, 2, 0 },
  { 0, 0, 0, 0, 0, 1, 0, 0 },
  { 0, 6, 0, 0, 5, 0, 0, 0 } };

#define MAXALPHA 20       /* maximal length of alphabet */

static short  alias[MAXALPHA + 1];
static int    pair[MAXALPHA + 1][MAXALPHA + 1];
/* rtype[pair[i][j]]:=pair[j][i] */
static int    rtype[8] = {
  0, 2, 1, 4, 3, 6, 5, 7
};

#ifdef _OPENMP
#pragma omp threadprivate(Law_and_Order, BP_pair, alias, pair, rtype)
#endif

/* for backward compatibility */
#define ENCODE(c) encode_char(c)

static INLINE int
encode_char(char c)
{
  /* return numerical representation of base used e.g. in pair[][] */
  int code;

  c = toupper(c);

  if (energy_set > 0) {
    code = (int)(c - 'A') + 1;
  } else {
    const char *pos;
    pos = strchr(Law_and_Order, c);
    if (pos == NULL)
      code = 0;
    else
      code = (int)(pos - Law_and_Order);

    if (code > 5)
      code = 0;

    if (code > 4)
      code--;           /* make T and U equivalent */
  }

  return code;
}


/*@+boolint +charint@*/
/*@null@*/
extern char *nonstandards;

static INLINE void
make_pair_matrix(void)
{
  int i, j;

  if (energy_set == 0) {
    for (i = 0; i < 5; i++)
      alias[i] = (short)i;
    alias[5]  = 3;  /* X <-> G */
    alias[6]  = 2;  /* K <-> C */
    alias[7]  = 0;  /* I <-> default base '@' */
    for (i = 0; i < NBASES; i++)
      for (j = 0; j < NBASES; j++)
        pair[i][j] = BP_pair[i][j];
    if (noGU)
      pair[3][4] = pair[4][3] = 0;

    if (nonstandards != NULL) {
      /* allow nonstandard bp's */
      for (i = 0; i < (int)strlen(nonstandards); i += 2)
        pair[encode_char(nonstandards[i])]
        [encode_char(nonstandards[i + 1])] = 7;
    }

    for (i = 0; i < NBASES; i++)
      for (j = 0; j < NBASES; j++)
        rtype[pair[i][j]] = pair[j][i];
  } else {
    for (i = 0; i <= MAXALPHA; i++)
      for (j = 0; j <= MAXALPHA; j++)
        pair[i][j] = 0;
    if (energy_set == 1) {
      for (i = 1; i < MAXALPHA; ) {
        alias[i++]  = 3;      /* A <-> G */
        alias[i++]  = 2;      /* B <-> C */
      }
      for (i = 1; i < MAXALPHA; i++) {
        pair[i][i + 1] = 2;       /* AB <-> GC */
        i++;
        pair[i][i - 1] = 1;       /* BA <-> CG */
      }
    } else if (energy_set == 2) {
      for (i = 1; i < MAXALPHA; ) {
        alias[i++]  = 1;      /* A <-> A*/
        alias[i++]  = 4;      /* B <-> U */
      }
      for (i = 1; i < MAXALPHA; i++) {
        pair[i][i + 1] = 5;       /* AB <-> AU */
        i++;
        pair[i][i - 1] = 6;       /* BA <-> UA */
      }
    } else if (energy_set == 3) {
      for (i = 1; i < MAXALPHA - 2; ) {
        alias[i++]  = 3;    /* A <-> G */
        alias[i++]  = 2;    /* B <-> C */
        alias[i++]  = 1;    /* C <-> A */
        alias[i++]  = 4;    /* D <-> U */
      }
      for (i = 1; i < MAXALPHA - 2; i++) {
        pair[i][i + 1] = 2;     /* AB <-> GC */
        i++;
        pair[i][i - 1] = 1;     /* BA <-> CG */
        i++;
        pair[i][i + 1] = 5;     /* CD <-> AU */
        i++;
        pair[i][i - 1] = 6;     /* DC <-> UA */
      }
    } else {
      vrna_message_error("What energy_set are YOU using??");
    }

    for (i = 0; i <= MAXALPHA; i++)
      for (j = 0; j <= MAXALPHA; j++)
        rtype[pair[i][j]] = pair[j][i];
  }
}


static INLINE short *
encode_sequence(const char  *sequence,
                short       how)
{
  unsigned int  i, l = (unsigned int)strlen(sequence);
  short         *S = (short *)vrna_alloc(sizeof(short) * (l + 2));

  switch (how) {
    /* standard encoding as always used for S */
    case 0:
      for (i = 1; i <= l; i++)    /* make numerical encoding of sequence */
        S[i] = (short)encode_char(sequence[i - 1]);
      S[l + 1]  = S[1];
      S[0]      = (short)l;
      break;
    /* encoding for mismatches of nostandard bases (normally used for S1) */
    case 1:
      for (i = 1; i <= l; i++)
        S[i] = alias[(short)encode_char(sequence[i - 1])];
      S[l + 1]  = S[1];
      S[0]      = S[l];
      break;
  }

  return S;
}


#endif /* VIENNA_RNA_PACKAGE_PAIR_MAT_H */