/*
  gquad.c

  Ronny Lorenz 2012

  Vienna RNA package
*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>

#include "config.h" // chzs
#include "fold_vars.h"
#include "data_structures.h"
#include "energy_const.h"
#include "utils.h"
#include "aln_util.h"
#include "gquad.h"

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

/**
 *  Use this macro to loop over each G-quadruplex
 *  delimited by a and b within the subsequence [c,d]
 */
#define FOR_EACH_GQUAD(a, b, c, d)  \
          for((a) = (d) - VRNA_GQUAD_MIN_BOX_SIZE + 1; (a) >= (c); (a)--)\
            for((b) = (a) + VRNA_GQUAD_MIN_BOX_SIZE - 1;\
                (b) <= MIN2((d), (a) + VRNA_GQUAD_MAX_BOX_SIZE - 1);\
                (b)++)

/**
 *  This macro does almost the same as FOR_EACH_GQUAD() but keeps
 *  the 5' delimiter fixed. 'b' is the 3' delimiter of the gquad,
 *  for gquads within subsequence [a,c] that have 5' delimiter 'a'
 */
#define FOR_EACH_GQUAD_AT(a, b, c)  \
          for((b) = (a) + VRNA_GQUAD_MIN_BOX_SIZE - 1;\
              (b) <= MIN2((c), (a) + VRNA_GQUAD_MAX_BOX_SIZE - 1);\
              (b)++)


/*
#################################
# PRIVATE FUNCTION DECLARATIONS #
#################################
*/

PRIVATE INLINE
int *
get_g_islands(short *S);

PRIVATE INLINE
int *
get_g_islands_sub(short *S, int i, int j);

/**
 *  IMPORTANT:
 *  If you don't know how to use this function, DONT'T USE IT!
 *
 *  The function pointer this function takes as argument is
 *  used for individual calculations with each g-quadruplex
 *  delimited by [i,j].
 *  The function it points to always receives as first 3 arguments
 *  position i, the stack size L and an array l[3] containing the
 *  individual linker sizes.
 *  The remaining 4 (void *) pointers of the callback function receive
 *  the parameters 'data', 'P', 'aux1' and 'aux2' and thus may be
 *  used to pass whatever data you like to.
 *  As the names of those parameters suggest the convention is that
 *  'data' should be used as a pointer where data is stored into,
 *  e.g the MFE or PF and the 'P' parameter should actually be a
 *  'paramT *' or 'pf_paramT *' type.
 *  However, what you actually pass obviously depends on the
 *  function the pointer is pointing to.
 *
 *  Although all of this may look like an overkill, it is found
 *  to be almost as fast as implementing g-quadruplex enumeration
 *  in each individual scenario, i.e. code duplication.
 *  Using this function, however, ensures that all g-quadruplex
 *  enumerations are absolutely identical.
 */
PRIVATE
void
process_gquad_enumeration(int *gg,
                          int i,
                          int j,
                          void (*f)(int, int, int *,
                                    void *, void *, void *, void *),
                          void *data,
                          void *P,
                          void *aux1,
                          void *aux2);

/**
 *  MFE callback for process_gquad_enumeration()
 */
PRIVATE
void
gquad_mfe(int i,
          int L,
          int *l,
          void *data,
          void *P,
          void *NA,
          void *NA2);

PRIVATE
void
gquad_mfe_pos(int i,
              int L,
              int *l,
              void *data,
              void *P,
              void *Lmfe,
              void *lmfe);

PRIVATE
void
gquad_pos_exhaustive( int i,
                      int L,
                      int *l,
                      void *data,
                      void *P,
                      void *Lex,
                      void *lex);

/**
 * Partition function callback for process_gquad_enumeration()
 */
PRIVATE
void
gquad_pf( int i,
          int L,
          int *l,
          void *data,
          void *P,
          void *NA,
          void *NA2);

/**
 * Partition function callback for process_gquad_enumeration()
 * in contrast to gquad_pf() it stores the stack size L and
 * the linker lengths l[3] of the g-quadruplex that dominates
 * the interval [i,j]
 * (FLT_OR_DBL *)data must be 0. on entry
 */
PRIVATE
void
gquad_pf_pos( int i,
              int L,
              int *l,
              void *data,
              void *pf,
              void *Lmax,
              void *lmax);

/**
 * MFE (alifold) callback for process_gquad_enumeration()
 */
PRIVATE
void
gquad_mfe_ali(int i,
              int L,
              int *l,
              void *data,
              void *P,
              void *S,
              void *n_seq);

/**
 * MFE (alifold) callback for process_gquad_enumeration()
 * with seperation of free energy and penalty contribution
 */
PRIVATE
void
gquad_mfe_ali_en( int i,
                  int L,
                  int *l,
                  void *data,
                  void *P,
                  void *S,
                  void *n_seq);

PRIVATE
void
gquad_interact( int i,
                int L,
                int *l,
                void *data,
                void *pf,
                void *index,
                void *NA2);

PRIVATE
void
gquad_count(int i,
            int L,
            int *l,
            void *data,
            void *NA,
            void *NA2,
            void *NA3);

PRIVATE
void
gquad_count_layers( int i,
                    int L,
                    int *l,
                    void *data,
                    void *NA,
                    void *NA2,
                    void *NA3);

/* other useful static functions */

PRIVATE
int
gquad_ali_penalty(int i,
                  int L,
                  int l[3],
                  const short **S,
                  paramT *P);

/*
#########################################
# BEGIN OF PUBLIC FUNCTION DEFINITIONS  #
#      (all available in RNAlib)        #
#########################################
*/

/********************************
  Here are the G-quadruplex energy
  contribution functions
*********************************/

PUBLIC int E_gquad( int L,
                    int l[3],
                    paramT *P){

  int i, c = INF;

  for(i=0;i<3;i++){
    if(l[i] > VRNA_GQUAD_MAX_LINKER_LENGTH) return c;
    if(l[i] < VRNA_GQUAD_MIN_LINKER_LENGTH) return c;
  }
  if(L > VRNA_GQUAD_MAX_STACK_SIZE) return c;
  if(L < VRNA_GQUAD_MIN_STACK_SIZE) return c;
  
  gquad_mfe(0, L, l,
            (void *)(&c),
            (void *)P,
            NULL,
            NULL);
  return c;
}

PUBLIC FLT_OR_DBL exp_E_gquad(int L,
                              int l[3],
                              pf_paramT *pf){

  int i;
  FLT_OR_DBL q = 0.;

  for(i=0;i<3;i++){
    if(l[i] > VRNA_GQUAD_MAX_LINKER_LENGTH) return q;
    if(l[i] < VRNA_GQUAD_MIN_LINKER_LENGTH) return q;
  }
  if(L > VRNA_GQUAD_MAX_STACK_SIZE) return q;
  if(L < VRNA_GQUAD_MIN_STACK_SIZE) return q;

  gquad_pf( 0, L, l,
            (void *)(&q),
            (void *)pf,
            NULL,
            NULL);
  return q;
}

PUBLIC int E_gquad_ali( int i,
                        int L,
                        int l[3],
                        const short **S,
                        int n_seq,
                        paramT *P){

  int en[2];
  E_gquad_ali_en(i, L, l, S, n_seq, en, P);
  return en[0] + en[1];
}


PUBLIC void E_gquad_ali_en( int i,
                            int L,
                            int l[3],
                            const short **S,
                            int n_seq,
                            int en[2],
                            paramT *P){

  int j;
  en[0] = en[1] = INF;

  for(j=0;j<3;j++){
    if(l[j] > VRNA_GQUAD_MAX_LINKER_LENGTH) return;
    if(l[j] < VRNA_GQUAD_MIN_LINKER_LENGTH) return;
  }
  if(L > VRNA_GQUAD_MAX_STACK_SIZE) return;
  if(L < VRNA_GQUAD_MIN_STACK_SIZE) return;

  gquad_mfe_ali_en( i, L, l,
                    (void *)(&(en[0])),
                    (void *)P,
                    (void *)S,
                    (void *)(&n_seq));
}

/********************************
  Now, the triangular matrix
  generators for the G-quadruplex
  contributions are following
*********************************/

PUBLIC int *get_gquad_matrix(short *S, paramT *P){

  int n, size, i, j, *gg, *my_index, *data;

  n         = S[0];
  my_index  = get_indx(n);
  gg        = get_g_islands(S);
  size      = (n * (n+1))/2 + 2;
  data      = (int *)space(sizeof(int) * size);

  /* prefill the upper triangular matrix with INF */
  for(i = 0; i < size; i++) data[i] = INF;

  FOR_EACH_GQUAD(i, j, 1, n){
    process_gquad_enumeration(gg, i, j,
                              &gquad_mfe,
                              (void *)(&(data[my_index[j]+i])),
                              (void *)P,
                              NULL,
                              NULL);
  }

  free(my_index);
  free(gg);
  return data;
}

PUBLIC FLT_OR_DBL *get_gquad_pf_matrix( short *S,
                                        FLT_OR_DBL *scale,
                                        pf_paramT *pf){

  int n, size, *gg, i, j, *my_index;
  FLT_OR_DBL *data;


  n         = S[0];
  size      = (n * (n+1))/2 + 2;
  data      = (FLT_OR_DBL *)space(sizeof(FLT_OR_DBL) * size);
  gg        = get_g_islands(S);
  my_index  = get_iindx(n);

  FOR_EACH_GQUAD(i, j, 1, n){
    process_gquad_enumeration(gg, i, j,
                              &gquad_pf,
                              (void *)(&(data[my_index[i]-j])),
                              (void *)pf,
                              NULL,
                              NULL);
    data[my_index[i]-j] *= scale[j-i+1];
  }

  free(my_index);
  free(gg);
  return data;
}

PUBLIC int *get_gquad_ali_matrix( short *S_cons,
                                  short **S,
                                  int n_seq,
                                  paramT *P){

  int n, size, *data, *gg;
  int i, j, *my_index;


  n         = S[0][0];
  size      = (n * (n+1))/2 + 2;
  data      = (int *)space(sizeof(int) * size);
  gg        = get_g_islands(S_cons);
  my_index  = get_indx(n);

  /* prefill the upper triangular matrix with INF */
  for(i=0;i<size;i++) data[i] = INF;

  FOR_EACH_GQUAD(i, j, 1, n){
    process_gquad_enumeration(gg, i, j,
                              &gquad_mfe_ali,
                              (void *)(&(data[my_index[j]+i])),
                              (void *)P,
                              (void *)S,
                              (void *)(&n_seq));
  }

  free(my_index);
  free(gg);
  return data;
}

PUBLIC int **get_gquad_L_matrix(short *S,
                                int start,
                                int maxdist,
                                int **g,
                                paramT *P){

  int **data;
  int n, i, j, k, l, *gg;
  
  n   = S[0];
  gg  = get_g_islands_sub(S, start, MIN2(n, start + maxdist + 4));

  if(g){ /* we just update the gquadruplex contribution for the current
            start and rotate the rest */
    data = g;
    /* we re-use the memory allocated previously */
    data[start] = data[start + maxdist + 5];
    data[start + maxdist + 5] = NULL;

    /* prefill with INF */
    for(i = 0; i < maxdist + 5; i++)
      data[start][i] = INF;

    /*  now we compute contributions for all gquads with 5' delimiter at
        position 'start'
    */
    FOR_EACH_GQUAD_AT(start, j, start + maxdist + 4){
      process_gquad_enumeration(gg, start, j,
                                &gquad_mfe,
                                (void *)(&(data[start][j-start])),
                                (void *)P,
                                NULL,
                                NULL);
    }

  } else { /* create a new matrix from scratch since this is the first
              call to this function */

    /* allocate memory and prefill with INF */
    data = (int **) space(sizeof(int *) * (n+1));
    for(k = n; (k>n-maxdist-5) && (k>=0); k--){
      data[k] = (int *) space(sizeof(int)*(maxdist+5));
      for(i = 0; i < maxdist+5; i++) data[k][i] = INF;
    }
    
    /* compute all contributions for the gquads in this interval */
    FOR_EACH_GQUAD(i, j, n - maxdist - 4, n){
      process_gquad_enumeration(gg, i, j,
                                &gquad_mfe,
                                (void *)(&(data[i][j-i])),
                                (void *)P,
                                NULL,
                                NULL);
    }
  }

  gg += start - 1;
  free(gg);
  return data;
}

PUBLIC plist *get_plist_gquad_from_db(const char *structure, float pr){
  int x, size, actual_size, L, n, ge, ee, gb, l[3];
  plist *pl;

  actual_size = 0;
  ge          = 0;
  n           = 2;
  size        = strlen(structure);
  pl          = (plist *)space(n*size*sizeof(plist));

  while((ee = parse_gquad(structure + ge, &L, l)) > 0){
    ge += ee;
    gb = ge - L*4 - l[0] - l[1] - l[2] + 1;
    /* add pseudo-base pair encloding gquad */
    for(x = 0; x < L; x++){
      if (actual_size >= n * size - 5){
        n *= 2;
        pl = (plist *)xrealloc(pl, n * size * sizeof(plist));
      }
      pl[actual_size].i = gb + x;
      pl[actual_size].j = ge + x - L + 1;
      pl[actual_size].p = pr;
      pl[actual_size++].type = 0;

      pl[actual_size].i = gb + x;
      pl[actual_size].j = gb + x + l[0] + L;
      pl[actual_size].p = pr;
      pl[actual_size++].type = 0;

      pl[actual_size].i = gb + x + l[0] + L;
      pl[actual_size].j = ge + x - 2*L - l[2] + 1;
      pl[actual_size].p = pr;
      pl[actual_size++].type = 0;

      pl[actual_size].i = ge + x - 2*L - l[2] + 1;
      pl[actual_size].j = ge + x - L + 1;
      pl[actual_size].p = pr;
      pl[actual_size++].type = 0;
    }
  } 

  pl[actual_size].i = pl[actual_size].j = 0;
  pl[actual_size++].p = 0;
  pl = (plist *)xrealloc(pl, actual_size * sizeof(plist));
  return pl;
}

PUBLIC void get_gquad_pattern_mfe(short *S,
                                  int i,
                                  int j,
                                  paramT *P,
                                  int *L,
                                  int l[3]){

  int *gg = get_g_islands_sub(S, i, j);
  int c = INF;

  process_gquad_enumeration(gg, i, j,
                            &gquad_mfe_pos,
                            (void *)(&c),
                            (void *)P,
                            (void *)L,
                            (void *)l);

  gg += i - 1;
  free(gg);
}

PUBLIC void
get_gquad_pattern_exhaustive( short *S,
                              int i,
                              int j,
                              paramT *P,
                              int *L,
                              int *l,
                              int threshold){

  int *gg = get_g_islands_sub(S, i, j);

  process_gquad_enumeration(gg, i, j,
                            &gquad_pos_exhaustive,
                            (void *)(&threshold),
                            (void *)P,
                            (void *)L,
                            (void *)l);

  gg += i - 1;
  free(gg);
}

PUBLIC void get_gquad_pattern_pf( short *S,
                                  int i,
                                  int j,
                                  pf_paramT *pf,
                                  int *L,
                                  int l[3]){

  int *gg = get_g_islands_sub(S, i, j);
  FLT_OR_DBL q = 0.;

  process_gquad_enumeration(gg, i, j,
                            &gquad_pf_pos,
                            (void *)(&q),
                            (void *)pf,
                            (void *)L,
                            (void *)l);

  gg += i - 1;
  free(gg);
}

PUBLIC plist *get_plist_gquad_from_pr(short *S,
                                      int gi,
                                      int gj,
                                      FLT_OR_DBL *G,
                                      FLT_OR_DBL *probs,
                                      FLT_OR_DBL *scale,
                                      pf_paramT *pf){

  int L, l[3];
  return  get_plist_gquad_from_pr_max(S, gi, gj, G, probs, scale, &L, l, pf);
}


PUBLIC plist *get_plist_gquad_from_pr_max(short *S,
                                      int gi,
                                      int gj,
                                      FLT_OR_DBL *G,
                                      FLT_OR_DBL *probs,
                                      FLT_OR_DBL *scale,
                                      int *Lmax,
                                      int lmax[3],
                                      pf_paramT *pf){ 

  int n, size, *gg, counter, i, j, *my_index;
  FLT_OR_DBL pp, *tempprobs;
  plist *pl;
  
  n         = S[0];
  size      = (n * (n + 1))/2 + 2;
  tempprobs = (FLT_OR_DBL *)space(sizeof(FLT_OR_DBL) * size);
  pl        = (plist *)space((S[0]*S[0])*sizeof(plist));
  gg        = get_g_islands_sub(S, gi, gj);
  counter   = 0;
  my_index  = get_iindx(n);

  process_gquad_enumeration(gg, gi, gj,
                            &gquad_interact,
                            (void *)tempprobs,
                            (void *)pf,
                            (void *)my_index,
                            NULL);

  pp = 0.;
  process_gquad_enumeration(gg, gi, gj,
                            &gquad_pf_pos,
                            (void *)(&pp),
                            (void *)pf,
                            (void *)Lmax,
                            (void *)lmax);

  pp = probs[my_index[gi]-gj] * scale[gj-gi+1] / G[my_index[gi]-gj];
  for (i=gi;i<gj; i++) {
    for (j=i; j<=gj; j++) {
      if (tempprobs[my_index[i]-j]>0.) {
        pl[counter].i=i;
        pl[counter].j=j;
        pl[counter++].p = pp * tempprobs[my_index[i]-j];
      }
    }
  }
  pl[counter].i = pl[counter].j = 0;
  pl[counter++].p = 0.;
  /* shrink memory to actual size needed */
  pl = (plist *) xrealloc(pl, counter * sizeof(plist));

  gg += gi - 1; free(gg);
  free(my_index);
  free (tempprobs);
  return pl;
}

PUBLIC int
get_gquad_count(short *S,
                int i,
                int j){

  int *gg     = get_g_islands_sub(S, i, j);
  int p,q,counter = 0;

  FOR_EACH_GQUAD(p, q, i, j)
    process_gquad_enumeration(gg, p, q,
                              &gquad_count,
                              (void *)(&counter),
                              NULL,
                              NULL,
                              NULL);

  gg += i - 1;
  free(gg);
  return counter;
}

PUBLIC int
get_gquad_layer_count(short *S,
                      int i,
                      int j){

  int *gg     = get_g_islands_sub(S, i, j);
  int p,q,counter = 0;

  FOR_EACH_GQUAD(p, q, i, j)
    process_gquad_enumeration(gg, p, q,
                              &gquad_count_layers,
                              (void *)(&counter),
                              NULL,
                              NULL,
                              NULL);

  gg += i - 1;
  free(gg);
  return counter;
}

PUBLIC int parse_gquad(const char *struc, int *L, int l[3]) {
  int i, il, start, end, len;

  for (i=0; struc[i] && struc[i]!='+'; i++);
  if (struc[i] == '+') { /* start of gquad */
    for (il=0; il<=3; il++) {
      start=i; /* pos of first '+' */
      while (struc[++i] == '+'){
        if((il) && (i-start == *L))
          break;
      }
      end=i; len=end-start; 
      if (il==0) *L=len;
      else if (len!=*L)
        nrerror("unequal stack lengths in gquad");
      if (il==3) break;
      while (struc[++i] == '.'); /* linker */
      l[il] = i-end;
      if (struc[i] != '+')
        nrerror("illegal character in gquad linker region");
    }
  }
  else return 0;
  /* printf("gquad at %d %d %d %d %d\n", end, *L, l[0], l[1], l[2]); */
  return end;
}



/*
#########################################
# BEGIN OF PRIVATE FUNCTION DEFINITIONS #
#          (internal use only)          #
#########################################
*/

PRIVATE int gquad_ali_penalty(int i,
                              int L,
                              int l[3],
                              const short **S,
                              paramT *P){

  int s, cnt;
  int penalty     = 0;
  int gg_mismatch = 0;

  /* check for compatibility in the alignment */
  for(s = 0; S[s]; s++){
    unsigned int  ld  = 0; /* !=0 if layer destruction was detected */
    int           pen = 0;

    /* check bottom layer */
    if(S[s][i] != 3)                            ld |= 1U;
    if(S[s][i + L + l[0]] != 3)                 ld |= 2U;
    if(S[s][i + 2*L + l[0] + l[1]] != 3)        ld |= 4U;
    if(S[s][i + 3*L + l[0] + l[1] + l[2]] != 3) ld |= 8U;
     /* add 1x penalty for missing bottom layer */
    if(ld) pen += VRNA_GQUAD_MISMATCH_PENALTY;

    /* check top layer */
    ld = 0;
    if(S[s][i + L - 1] != 3)                        ld |= 1U;
    if(S[s][i + 2*L + l[0] - 1] != 3)               ld |= 2U;
    if(S[s][i + 3*L + l[0] + l[1] - 1] != 3)        ld |= 4U;
    if(S[s][i + 4*L + l[0] + l[1] + l[2] - 1] != 3) ld |= 8U;
     /* add 1x penalty for missing top layer */
    if(ld) pen += VRNA_GQUAD_MISMATCH_PENALTY;

    /* check inner layers */
    for(cnt=1;cnt<L-1;cnt++){
      if(S[s][i + cnt] != 3)                            ld |= 1U;
      if(S[s][i + L + l[0] + cnt] != 3)                 ld |= 2U;
      if(S[s][i + 2*L + l[0] + l[1] + cnt] != 3)        ld |= 4U;
      if(S[s][i + 3*L + l[0] + l[1] + l[2] + cnt] != 3) ld |= 8U;
      /* add 2x penalty for missing inner layer */
      if(ld) pen += 2*VRNA_GQUAD_MISMATCH_PENALTY;
    }

    /* if all layers are missing, we have a complete gg mismatch */
    if(pen >= (2*VRNA_GQUAD_MISMATCH_PENALTY * (L-1)))
      gg_mismatch++;

    /* add the penalty to the score */
    penalty += pen;
  }
  /* if gg_mismatch exceeds maximum allowed, this g-quadruplex is forbidden */
  if(gg_mismatch > VRNA_GQUAD_MISMATCH_NUM_ALI) return INF;
  else return penalty;
}


PRIVATE void gquad_mfe( int i,
                        int L,
                        int *l,
                        void *data,
                        void *P,
                        void *NA,
                        void *NA2){

  int cc = ((paramT *)P)->gquad[L][l[0] + l[1] + l[2]];
  if(cc < *((int *)data))
    *((int *)data) = cc;
}

PRIVATE void gquad_mfe_pos( int i,
                            int L,
                            int *l,
                            void *data,
                            void *P,
                            void *Lmfe,
                            void *lmfe){

  int cc = ((paramT *)P)->gquad[L][l[0] + l[1] + l[2]];
  if(cc < *((int *)data)){
    *((int *)data)        = cc;
    *((int *)Lmfe)        = L;
    *((int *)lmfe)        = l[0];
    *(((int *)lmfe) + 1)  = l[1];
    *(((int *)lmfe) + 2)  = l[2];
  }
}

PRIVATE
void
gquad_pos_exhaustive( int i,
                      int L,
                      int *l,
                      void *data,
                      void *P,
                      void *Lex,
                      void *lex){

  int cnt;
  int cc = ((paramT *)P)->gquad[L][l[0] + l[1] + l[2]];
  if(cc <= *((int *)data)){
    /*  since Lex is an array of L values and lex an
        array of l triples we need to find out where
        the current gquad position is to be stored...
		the below implementation might be slow but we
		still use it for now
    */
    for(cnt = 0; ((int *)Lex)[cnt] != -1; cnt++);

    *((int *)Lex + cnt)           = L;
    *((int *)Lex + cnt + 1)       = -1;
    *(((int *)lex) + (3*cnt) + 0) = l[0];
    *(((int *)lex) + (3*cnt) + 1) = l[1];
    *(((int *)lex) + (3*cnt) + 2) = l[2];
  }
}

PRIVATE
void
gquad_count(int i,
            int L,
            int *l,
            void *data,
            void *NA,
            void *NA2,
            void *NA3){

  *((int *)data) += 1;
}

PRIVATE
void
gquad_count_layers( int i,
                    int L,
                    int *l,
                    void *data,
                    void *NA,
                    void *NA2,
                    void *NA3){

  *((int *)data) += L;
}


PRIVATE void gquad_pf(int i,
                      int L,
                      int *l,
                      void *data,
                      void *pf,
                      void *NA,
                      void *NA2){

  *((FLT_OR_DBL *)data) += ((pf_paramT *)pf)->expgquad[L][l[0] + l[1] + l[2]];
}

PRIVATE void gquad_pf_pos(int i,
                          int L,
                          int *l,
                          void *data,
                          void *pf,
                          void *Lmax,
                          void *lmax){

  FLT_OR_DBL gq = ((pf_paramT *)pf)->expgquad[L][l[0] + l[1] + l[2]];
  if(gq > *((FLT_OR_DBL *)data)){
    *((FLT_OR_DBL *)data) = gq;
    *((int *)Lmax)        = L;
    *((int *)lmax)        = l[0];
    *(((int *)lmax) + 1)  = l[1];
    *(((int *)lmax) + 2)  = l[2];
  }
}

PRIVATE void gquad_mfe_ali( int i,
                            int L,
                            int *l,
                            void *data,
                            void *P,
                            void *S,
                            void *n_seq){

  int en[2], cc;
  en[0] = en[1] = INF;
  gquad_mfe_ali_en(i, L, l, (void *)(&(en[0])), P, S, n_seq);
  if(en[1] != INF){
    cc  = en[0] + en[1];
    if(cc < *((int *)data)) *((int *)data) = cc;
  }
}

PRIVATE void gquad_mfe_ali_en(int i,
                              int L,
                              int *l,
                              void *data,
                              void *P,
                              void *S,
                              void *n_seq){

  int en[2], cc, dd;
  en[0] = ((paramT *)P)->gquad[L][l[0] + l[1] + l[2]] * (*(int *)n_seq);
  en[1] = gquad_ali_penalty(i, L, l, (const short **)S, (paramT *)P);
  if(en[1] != INF){
    cc = en[0] + en[1];
    dd = ((int *)data)[0] + ((int *)data)[1];
    if(cc < dd){
      ((int *)data)[0] = en[0];
      ((int *)data)[1] = en[1];
    }
  }
}

PRIVATE void gquad_interact(int i,
                      int L,
                      int *l,
                      void *data,
                      void *pf,
                      void *index,
                      void *NA2){

  int x, *idx;
  FLT_OR_DBL gq, *pp;

  idx = (int *)index;
  pp  = (FLT_OR_DBL *)data;
  gq  = exp_E_gquad(L, l, (pf_paramT *)pf);

  for(x = 0; x < L; x++){
    pp[idx[i + x] - (i + x + 3*L + l[0] + l[1] + l[2])] += gq;
    pp[idx[i + x] - (i + x + L + l[0])] += gq;
    pp[idx[i + x + L + l[0]] - (i + x + 2*L + l[0] + l[1])] += gq;
    pp[idx[i + x + 2*L + l[0] + l[1]] - (i + x + 3*L + l[0] + l[1] + l[2])] += gq;
  }
  
}

PRIVATE INLINE int *get_g_islands(short *S){
  return get_g_islands_sub(S, 1, S[0]);
}

PRIVATE INLINE int *get_g_islands_sub(short *S, int i, int j){
  int x, *gg;

  gg = (int *)space(sizeof(int)*(j-i+2));
  gg -= i - 1;

  if(S[j]==3) gg[j] = 1;
  for(x = j - 1; x >= i; x--)
    if(S[x] == 3)
      gg[x] = gg[x+1]+1;

  return gg;
}

/**
 *  We could've also created a macro that loops over all G-quadruplexes
 *  delimited by i and j. However, for the fun of it we use this function
 *  that receives a pointer to a callback function which in turn does the
 *  actual computation for each quadruplex found.
 */
PRIVATE
void
process_gquad_enumeration(int *gg,
                          int i,
                          int j,
                          void (*f)(int, int, int *,
                                    void *, void *, void *, void *),
                          void *data,
                          void *P,
                          void *aux1,
                          void *aux2){

  int L, l[3], n, max_linker, maxl0, maxl1;

  n = j - i + 1;

  if((n >= VRNA_GQUAD_MIN_BOX_SIZE) && (n <= VRNA_GQUAD_MAX_BOX_SIZE))
    for(L = MIN2(gg[i], VRNA_GQUAD_MAX_STACK_SIZE);
        L >= VRNA_GQUAD_MIN_STACK_SIZE;
        L--)
      if(gg[j-L+1] >= L){
        max_linker = n-4*L;
        if(     (max_linker >= 3*VRNA_GQUAD_MIN_LINKER_LENGTH)
            &&  (max_linker <= 3*VRNA_GQUAD_MAX_LINKER_LENGTH)){
          maxl0 = MIN2( VRNA_GQUAD_MAX_LINKER_LENGTH,
                        max_linker - 2*VRNA_GQUAD_MIN_LINKER_LENGTH
                      );
          for(l[0] = VRNA_GQUAD_MIN_LINKER_LENGTH;
              l[0] <= maxl0;
              l[0]++)
            if(gg[i+L+l[0]] >= L){
              maxl1 = MIN2( VRNA_GQUAD_MAX_LINKER_LENGTH,
                            max_linker - l[0] - VRNA_GQUAD_MIN_LINKER_LENGTH
                          );
              for(l[1] = VRNA_GQUAD_MIN_LINKER_LENGTH;
                  l[1] <= maxl1;
                  l[1]++)
                if(gg[i + 2*L + l[0] + l[1]] >= L){
                  l[2] = max_linker - l[0] - l[1];
                  f(i, L, &(l[0]), data, P, aux1, aux2);
                }
            }
        }
      }
}