/* Last changed Time-stamp: <2009-02-24 15:17:17 ivo> */
/*
                  minimum free energy folding
                  for a set of aligned sequences

                  c Ivo Hofacker

                  Vienna RNA package
*/

/**
*** \file alifold.c
**/


#include <config.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <ctype.h>
#include <string.h>
#include <limits.h>

#include "fold.h"
#include "utils.h"
#include "energy_par.h"
#include "fold_vars.h"
#include "pair_mat.h"
#include "params.h"
#include "ribo.h"
#include "aln_util.h"
#include "loop_energies.h"
#include "gquad.h"
#include "alifold.h"

#ifdef _OPENMP
#include <omp.h>
#endif

#define PAREN
#define STACK_BULGE1  1     /* stacking energies for bulges of size 1 */
#define NEW_NINIO     1     /* new asymetry penalty */
#define MAXSECTORS    500   /* dimension for a backtrack array */
#define LOCALITY        0.  /* locality parameter for base-pairs */
#define UNIT 100
#define MINPSCORE -2 * UNIT

/*
#################################
# GLOBAL VARIABLES              #
#################################
*/
PUBLIC double cv_fact=1.; /* should be made static to not interfere with other threads */
PUBLIC double nc_fact=1.; /* should be made static to not interfere with other threads */

/*
#################################
# PRIVATE VARIABLES             #
#################################
*/
PRIVATE short           **S     = NULL;
PRIVATE short           **S5    = NULL;     /*S5[s][i] holds next base 5' of i in sequence s*/
PRIVATE short           **S3    = NULL;     /*Sl[s][i] holds next base 3' of i in sequence s*/
PRIVATE char            **Ss    = NULL;
PRIVATE unsigned short  **a2s   = NULL;
PRIVATE paramT          *P      = NULL;
PRIVATE int             *indx   = NULL;     /* index for moving in the triangle matrices c[] and fMl[]*/
PRIVATE int             *c      = NULL;     /* energy array, given that i-j pair */
PRIVATE int             *cc     = NULL;     /* linear array for calculating canonical structures */
PRIVATE int             *cc1    = NULL;     /*   "     "        */
PRIVATE int             *f5     = NULL;     /* energy of 5' end */
PRIVATE int             *fML    = NULL;     /* multi-loop auxiliary energy array */
PRIVATE int             *Fmi    = NULL;     /* holds row i of fML (avoids jumps in memory) */
PRIVATE int             *DMLi   = NULL;     /* DMLi[j] holds MIN(fML[i,k]+fML[k+1,j])  */
PRIVATE int             *DMLi1  = NULL;     /*             MIN(fML[i+1,k]+fML[k+1,j])  */
PRIVATE int             *DMLi2  = NULL;     /*             MIN(fML[i+2,k]+fML[k+1,j])  */
PRIVATE int             *pscore = NULL;     /* precomputed array of pair types */
PRIVATE int             init_length = -1;
PRIVATE sect            sector[MAXSECTORS]; /* stack of partial structures for backtracking */
PRIVATE bondT           *base_pair2 = NULL;
PRIVATE int             circular    = 0;
PRIVATE int             with_gquad  = 0;

PRIVATE int             *ggg        = NULL; /* minimum free energies of the gquadruplexes */
PRIVATE char            *cons_seq   = NULL;
PRIVATE short           *S_cons     = NULL;

#ifdef _OPENMP

#pragma omp threadprivate(S, S5, S3, Ss, a2s, P, indx, c, cc, cc1, f5, fML, Fmi, DMLi, DMLi1, DMLi2,\
                          pscore, init_length, sector, base_pair2,\
                          ggg, with_gquad, cons_seq, S_cons)

#endif

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

PRIVATE void  init_alifold(int length);
PRIVATE void  get_arrays(unsigned int size);
PRIVATE void  make_pscores(const short *const *S, const char **AS, int n_seq, const char *structure);
PRIVATE int   fill_arrays(const char **strings);
PRIVATE void  backtrack(const char **strings, int s);

PRIVATE void    energy_of_alistruct_pt(const char **sequences,short * ptable, int n_seq, int *energy);
PRIVATE void    stack_energy_pt(int i, const char **sequences, short *ptable,  int n_seq, int *energy);
PRIVATE int     ML_Energy_pt(int i, int n_seq, short *pt);
PRIVATE int     EL_Energy_pt(int i, int n_seq, short *pt);

PRIVATE void    en_corr_of_loop_gquad(int i,
                                      int j,
                                      const char **sequences,
                                      const char *structure,
                                      short *pt,
                                      int *loop_idx,
                                      int n_seq,
                                      int en[2]);

/*
#################################
# BEGIN OF FUNCTION DEFINITIONS #
#################################
*/

/* unsafe function that will be replaced by a threadsafe companion in the future */
PRIVATE void init_alifold(int length){

#ifdef _OPENMP
/* Explicitly turn off dynamic threads */
  omp_set_dynamic(0);
#endif

  if (length < 1) nrerror("initialize_fold: argument must be greater 0");
  free_alifold_arrays();
  get_arrays((unsigned) length);
  init_length = length;

  indx = get_indx((unsigned)length);

  update_alifold_params();
}

PRIVATE void get_arrays(unsigned int size){
  if(size >= (unsigned int)sqrt((double)INT_MAX))
    nrerror("get_arrays@alifold.c: sequence length exceeds addressable range");

  c       = (int *) space(sizeof(int)*((size*(size+1))/2+2));
  fML     = (int *) space(sizeof(int)*((size*(size+1))/2+2));
  pscore  = (int *) space(sizeof(int)*((size*(size+1))/2+2));
  f5      = (int *) space(sizeof(int)*(size+2));
  cc      = (int *) space(sizeof(int)*(size+2));
  cc1     = (int *) space(sizeof(int)*(size+2));
  Fmi     = (int *) space(sizeof(int)*(size+1));
  DMLi    = (int *) space(sizeof(int)*(size+1));
  DMLi1   = (int *) space(sizeof(int)*(size+1));
  DMLi2   = (int *) space(sizeof(int)*(size+1));
  if(base_pair2) free(base_pair2);

  base_pair2 = (bondT *) space(sizeof(bondT)*(1+size/2+200)); /* add a guess of how many G's may be involved in a G quadruplex */
}

PUBLIC  void  free_alifold_arrays(void){
  if(indx)        free(indx);
  if(c)           free(c);
  if(fML)         free(fML);
  if(f5)          free(f5);
  if(cc)          free(cc);
  if(cc1)         free(cc1);
  if(pscore)      free(pscore);
  if(base_pair2)  free(base_pair2);
  if(Fmi)         free(Fmi);
  if(DMLi)        free(DMLi);
  if(DMLi1)       free(DMLi1);
  if(DMLi2)       free(DMLi2);
  if(P)           free(P);
  if(ggg)         free(ggg);
  if(cons_seq)    free(cons_seq);
  if(S_cons)      free(S_cons);

  indx = c = fML = f5 = cc = cc1 = Fmi = DMLi = DMLi1 = DMLi2 = ggg = NULL;
  pscore      = NULL;
  base_pair   = NULL;
  base_pair2  = NULL;
  P           = NULL;
  init_length = 0;
  cons_seq    = NULL;
  S_cons      = NULL; 
}


PUBLIC void alloc_sequence_arrays(const char **sequences, short ***S, short ***S5, short ***S3, unsigned short ***a2s, char ***Ss, int circ){
  unsigned int s, n_seq, length;
  if(sequences[0] != NULL){
    length = strlen(sequences[0]);
    for (s=0; sequences[s] != NULL; s++);
    n_seq = s;
    *S    = (short **)          space((n_seq+1) * sizeof(short *));
    *S5   = (short **)          space((n_seq+1) * sizeof(short *));
    *S3   = (short **)          space((n_seq+1) * sizeof(short *));
    *a2s  = (unsigned short **) space((n_seq+1) * sizeof(unsigned short *));
    *Ss   = (char **)           space((n_seq+1) * sizeof(char *));
    for (s=0; s<n_seq; s++) {
      if(strlen(sequences[s]) != length) nrerror("uneqal seqence lengths");
      (*S5)[s]  = (short *)         space((length + 2) * sizeof(short));
      (*S3)[s]  = (short *)         space((length + 2) * sizeof(short));
      (*a2s)[s] = (unsigned short *)space((length + 2) * sizeof(unsigned short));
      (*Ss)[s]  = (char *)          space((length + 2) * sizeof(char));
      (*S)[s]   = (short *)         space((length + 2) * sizeof(short));
      encode_ali_sequence(sequences[s], (*S)[s], (*S5)[s], (*S3)[s], (*Ss)[s], (*a2s)[s], circ);
    }
    (*S5)[n_seq]  = NULL;
    (*S3)[n_seq]  = NULL;
    (*a2s)[n_seq] = NULL;
    (*Ss)[n_seq]  = NULL;
    (*S)[n_seq]   = NULL;
  }
  else nrerror("alloc_sequence_arrays: no sequences in the alignment!");
}

PUBLIC void free_sequence_arrays(unsigned int n_seq, short ***S, short ***S5, short ***S3, unsigned short ***a2s, char ***Ss){
  unsigned int s;
  for (s=0; s<n_seq; s++) {
    free((*S)[s]);
    free((*S5)[s]);
    free((*S3)[s]);
    free((*a2s)[s]);
    free((*Ss)[s]);
  }
  free(*S);   *S    = NULL;
  free(*S5);  *S5   = NULL;
  free(*S3);  *S3   = NULL;
  free(*a2s); *a2s  = NULL;
  free(*Ss);  *Ss   = NULL;
}

PUBLIC void update_alifold_params(void){
  if(P) free(P);
  P = scale_parameters();
  make_pair_matrix();
  if (init_length < 0) init_length=0;
}

PUBLIC float alifold(const char **strings, char *structure){
  int  length, energy, s, n_seq;

  circular = 0;
  length = (int) strlen(strings[0]);

#ifdef _OPENMP
  /* always init everything since all global static variables are uninitialized when entering a thread */
  init_alifold(length);
#else
  if (length>init_length) init_alifold(length);
#endif
  if (fabs(P->temperature - temperature)>1e-6)  update_alifold_params();

  for (s=0; strings[s]!=NULL; s++);
  n_seq       = s;
  with_gquad  = P->model_details.gquad;

  alloc_sequence_arrays(strings, &S, &S5, &S3, &a2s, &Ss, circular);
  make_pscores((const short **) S, strings, n_seq, structure);

  energy = fill_arrays((const char **)strings);

  backtrack((const char **)strings, 0);

#ifdef PAREN
  parenthesis_structure(structure, base_pair2, length);
#else
  letter_structure(structure, base_pair2, length);
#endif

  /*
  *  Backward compatibility:
  *  This block may be removed if deprecated functions
  *  relying on the global variable "base_pair" vanishs from within the package!
  */
  base_pair = base_pair2;
  /*
  {
    if(base_pair) free(base_pair);
    base_pair = (bondT *)space(sizeof(bondT) * (1+length/2));
    memcpy(base_pair, base_pair2, sizeof(bondT) * (1+length/2));
  }
  */
  free_sequence_arrays(n_seq, &S, &S5, &S3, &a2s, &Ss);

  if (backtrack_type=='C')
    return (float) c[indx[length]+1]/(n_seq*100.);
  else if (backtrack_type=='M')
    return (float) fML[indx[length]+1]/(n_seq*100.);
  else
    return (float) f5[length]/(n_seq*100.);
}


/**
*** the actual forward recursion to fill the energy arrays
**/
PRIVATE int fill_arrays(const char **strings) {
  int   i, j, k, p, q, length, energy, new_c;
  int   decomp, MLenergy, new_fML;
  int   s, n_seq, *type, type_2, tt;

  /* count number of sequences */
  for (n_seq=0; strings[n_seq]!=NULL; n_seq++);

  type = (int *) space(n_seq*sizeof(int));
  length = strlen(strings[0]);

  /* init energies */

  if(with_gquad){
    cons_seq = consensus(strings);
    /* make g-island annotation of the consensus */
    S_cons = encode_sequence(cons_seq, 0);
    ggg = get_gquad_ali_matrix(S_cons, S, n_seq, P);
  }

  for (j=1; j<=length; j++){
    Fmi[j]=DMLi[j]=DMLi1[j]=DMLi2[j]=INF;
    for (i=(j>TURN?(j-TURN):1); i<j; i++) {
      c[indx[j]+i] = fML[indx[j]+i] = INF;
    }
  }

  /* begin recursions */
  for (i = length-TURN-1; i >= 1; i--) { /* i,j in [1..length] */
    for (j = i+TURN+1; j <= length; j++) {
      int ij, psc, l1, maxq, minq, up, c0;
      ij = indx[j]+i;

      for (s=0; s<n_seq; s++) {
        type[s] = pair[S[s][i]][S[s][j]];
        if (type[s]==0) type[s]=7;
      }

      psc = pscore[indx[j]+i];
      if (psc>=MINPSCORE) {   /* a pair to consider */
        int stackEnergy = INF;
        /* hairpin ----------------------------------------------*/


        for (new_c=s=0; s<n_seq; s++) {
          if ((a2s[s][j-1]-a2s[s][i])<3) new_c+=600;
          else  new_c += E_Hairpin(a2s[s][j-1]-a2s[s][i],type[s],S3[s][i],S5[s][j],Ss[s]+(a2s[s][i-1]), P);
       }
        /*--------------------------------------------------------
          check for elementary structures involving more than one
          closing pair.
          --------------------------------------------------------*/

        for (p = i+1; p <= MIN2(j-2-TURN,i+MAXLOOP+1) ; p++) {
          minq = j-i+p-MAXLOOP-2;
          if (minq<p+1+TURN) minq = p+1+TURN;
          for (q = minq; q < j; q++) {
            if (pscore[indx[q]+p]<MINPSCORE) continue;
            for (energy = s=0; s<n_seq; s++) {
              type_2 = pair[S[s][q]][S[s][p]]; /* q,p not p,q! */
              if (type_2 == 0) type_2 = 7;
              energy += E_IntLoop(a2s[s][p-1]-a2s[s][i], a2s[s][j-1]-a2s[s][q], type[s], type_2,
                                   S3[s][i], S5[s][j],
                                   S5[s][p], S3[s][q], P);
            }
            new_c = MIN2(new_c, energy + c[indx[q]+p]);
            if ((p==i+1)&&(j==q+1)) stackEnergy = energy; /* remember stack energy */
          } /* end q-loop */
        } /* end p-loop */

        /* multi-loop decomposition ------------------------*/
        decomp = DMLi1[j-1];
        if(dangles){
          for(s=0; s<n_seq; s++){
            tt = rtype[type[s]];
            decomp += E_MLstem(tt, S5[s][j], S3[s][i], P);
          }
        }
        else{
          for(s=0; s<n_seq; s++){
            tt = rtype[type[s]];
            decomp += E_MLstem(tt, -1, -1, P);
          }
        }
        MLenergy = decomp + n_seq*P->MLclosing;
        new_c = MIN2(new_c, MLenergy);

        if(with_gquad){
          decomp = 0;
          for(s=0;s<n_seq;s++){
            tt = type[s];
            if(dangles == 2)
              decomp += P->mismatchI[tt][S3[s][i]][S5[s][j]];
            if(tt > 2)
              decomp += P->TerminalAU;
          }
          for(p = i + 2; p < j - VRNA_GQUAD_MIN_BOX_SIZE; p++){
            l1    = p - i - 1;
            if(l1>MAXLOOP) break;
            if(S_cons[p] != 3) continue;
            minq  = j - i + p - MAXLOOP - 2;
            c0    = p + VRNA_GQUAD_MIN_BOX_SIZE - 1;
            minq  = MAX2(c0, minq);
            c0    = j - 1;
            maxq  = p + VRNA_GQUAD_MAX_BOX_SIZE + 1;
            maxq  = MIN2(c0, maxq);
            for(q = minq; q < maxq; q++){
              if(S_cons[q] != 3) continue;
              c0    = decomp + ggg[indx[q] + p] + n_seq * P->internal_loop[l1 + j - q - 1];
              new_c = MIN2(new_c, c0);
            }
          }

          p = i + 1;
          if(S_cons[p] == 3){
            if(p < j - VRNA_GQUAD_MIN_BOX_SIZE){
              minq  = j - i + p - MAXLOOP - 2;
              c0    = p + VRNA_GQUAD_MIN_BOX_SIZE - 1;
              minq  = MAX2(c0, minq);
              c0    = j - 3;
              maxq  = p + VRNA_GQUAD_MAX_BOX_SIZE + 1;
              maxq  = MIN2(c0, maxq);
              for(q = minq; q < maxq; q++){
                if(S_cons[q] != 3) continue;
                c0  = decomp + ggg[indx[q] + p] + n_seq * P->internal_loop[j - q - 1];
                new_c   = MIN2(new_c, c0);
              }
            }
          }
          q = j - 1;
          if(S_cons[q] == 3)
            for(p = i + 4; p < j - VRNA_GQUAD_MIN_BOX_SIZE; p++){
              l1    = p - i - 1;
              if(l1>MAXLOOP) break;
              if(S_cons[p] != 3) continue;
              c0  = decomp + ggg[indx[q] + p] + n_seq * P->internal_loop[l1];
              new_c   = MIN2(new_c, c0);
            }
        }

        new_c = MIN2(new_c, cc1[j-1]+stackEnergy);

        cc[j] = new_c - psc; /* add covariance bonnus/penalty */
        if (noLonelyPairs)
          c[ij] = cc1[j-1]+stackEnergy-psc;
        else
          c[ij] = cc[j];

      } /* end >> if (pair) << */

      else c[ij] = INF;
      /* done with c[i,j], now compute fML[i,j] */
      /* free ends ? -----------------------------------------*/

      new_fML = fML[ij+1]+n_seq*P->MLbase;
      new_fML = MIN2(fML[indx[j-1]+i]+n_seq*P->MLbase, new_fML);
      energy = c[ij];
      if(dangles){
        for (s=0; s<n_seq; s++) {
          energy += E_MLstem(type[s], S5[s][i], S3[s][j], P);
        }
      }
      else{
        for (s=0; s<n_seq; s++) {
          energy += E_MLstem(type[s], -1, -1, P);
        }
      }
      new_fML = MIN2(energy, new_fML);

      if(with_gquad){
        decomp = ggg[indx[j] + i] + n_seq * E_MLstem(0, -1, -1, P);
        new_fML = MIN2(new_fML, decomp);
      }

      /* modular decomposition -------------------------------*/
      for (decomp = INF, k = i+1+TURN; k <= j-2-TURN; k++)
        decomp = MIN2(decomp, Fmi[k]+fML[indx[j]+k+1]);

      DMLi[j] = decomp;               /* store for use in ML decompositon */
      new_fML = MIN2(new_fML,decomp);

      /* coaxial stacking deleted */

      fML[ij] = Fmi[j] = new_fML;     /* substring energy */
    } /* END for j */

    {
      int *FF; /* rotate the auxilliary arrays */
      FF = DMLi2; DMLi2 = DMLi1; DMLi1 = DMLi; DMLi = FF;
      FF = cc1; cc1=cc; cc=FF;
      for (j=1; j<=length; j++) {cc[j]=Fmi[j]=DMLi[j]=INF; }
    }
  } /* END for i */
  /* calculate energies of 5' and 3' fragments */

  f5[TURN + 1] = 0;
  switch(dangles){
    case 0:   for(j = TURN + 2; j <= length; j++){
                f5[j] = f5[j-1];
                if (c[indx[j]+1]<INF){
                  energy = c[indx[j]+1];
                  for(s = 0; s < n_seq; s++){
                    tt = pair[S[s][1]][S[s][j]];
                    if(tt==0) tt=7;
                    energy += E_ExtLoop(tt, -1, -1, P);
                  }
                  f5[j] = MIN2(f5[j], energy);
                }

                if(with_gquad){
                  if(ggg[indx[j]+1] < INF)
                    f5[j] = MIN2(f5[j], ggg[indx[j]+1]);
                }

                for(i = j - TURN - 1; i > 1; i--){
                  if(c[indx[j]+i]<INF){
                    energy = f5[i-1] + c[indx[j]+i];
                    for(s = 0; s < n_seq; s++){
                      tt = pair[S[s][i]][S[s][j]];
                      if(tt==0) tt=7;
                      energy += E_ExtLoop(tt, -1, -1, P);
                    }
                    f5[j] = MIN2(f5[j], energy);
                  }

                  if(with_gquad){
                    if(ggg[indx[j]+i] < INF)
                      f5[j] = MIN2(f5[j], f5[i-1] + ggg[indx[j]+i]);
                  }

                }
              }
              break;
    default:  for(j = TURN + 2; j <= length; j++){
                f5[j] = f5[j-1];
                if (c[indx[j]+1]<INF) {
                  energy = c[indx[j]+1];
                  for(s = 0; s < n_seq; s++){
                    tt = pair[S[s][1]][S[s][j]];
                    if(tt==0) tt=7;
                    energy += E_ExtLoop(tt, -1, (j<length) ? S3[s][j] : -1, P);
                  }
                  f5[j] = MIN2(f5[j], energy);
                }

                if(with_gquad){
                  if(ggg[indx[j]+1] < INF)
                    f5[j] = MIN2(f5[j], ggg[indx[j]+1]);
                }

                for(i = j - TURN - 1; i > 1; i--){
                  if (c[indx[j]+i]<INF) {
                    energy = f5[i-1] + c[indx[j]+i];
                    for(s = 0; s < n_seq; s++){
                      tt = pair[S[s][i]][S[s][j]];
                      if(tt==0) tt=7;
                      energy += E_ExtLoop(tt, S5[s][i], (j < length) ? S3[s][j] : -1, P);
                    }
                    f5[j] = MIN2(f5[j], energy);
                  }

                  if(with_gquad){
                    if(ggg[indx[j]+i] < INF)
                      f5[j] = MIN2(f5[j], f5[i-1] + ggg[indx[j]+i]);
                  }

                }
              }
              break;
  }
  free(type);
  return(f5[length]);
}

#include "alicircfold.inc"

/**
*** backtrack in the energy matrices to obtain a structure with MFE
**/
PRIVATE void backtrack(const char **strings, int s) {
  /*------------------------------------------------------------------
    trace back through the "c", "f5" and "fML" arrays to get the
    base pairing list. No search for equivalent structures is done.
    This inverts the folding procedure, hence it's very fast.
    ------------------------------------------------------------------*/
   /* normally s=0.
     If s>0 then s items have been already pushed onto the sector stack */
  int   i, j, k, p, q, length, energy, up, c0, l1, minq, maxq;
  int   type_2, tt, mm;
  int   b=0, cov_en = 0;
  int   n_seq;
  int *type;
  length = strlen(strings[0]);
  for (n_seq=0; strings[n_seq]!=NULL; n_seq++);
  type = (int *) space(n_seq*sizeof(int));

  if (s==0) {
    sector[++s].i = 1;
    sector[s].j = length;
    sector[s].ml = (backtrack_type=='M') ? 1 : ((backtrack_type=='C')?2:0);
  }
  while (s>0) {
    int ss, ml, fij, fi, cij, traced, i1, j1, d3, d5, jj=0, gq=0;
    int canonical = 1;     /* (i,j) closes a canonical structure */
    i  = sector[s].i;
    j  = sector[s].j;
    ml = sector[s--].ml;   /* ml is a flag indicating if backtracking is to
                              occur in the fML- (1) or in the f-array (0) */
    if (ml==2) {
      base_pair2[++b].i = i;
      base_pair2[b].j   = j;
      cov_en += pscore[indx[j]+i];
      goto repeat1;
    }

    if (j < i+TURN+1) continue; /* no more pairs in this interval */

    fij = (ml)? fML[indx[j]+i] : f5[j];
    fi  = (ml)?(fML[indx[j-1]+i]+n_seq*P->MLbase):f5[j-1];

    if (fij == fi) {  /* 3' end is unpaired */
      sector[++s].i = i;
      sector[s].j   = j-1;
      sector[s].ml  = ml;
      continue;
    }

    if (ml == 0) { /* backtrack in f5 */
      switch(dangles){
        case 0:   /* j or j-1 is paired. Find pairing partner */
                  for (i=j-TURN-1,traced=0; i>=1; i--) {
                    int cc, en;
                    jj = i-1;
                    if (c[indx[j]+i]<INF) {
                      en = c[indx[j]+i] + f5[i-1];
                      for(ss = 0; ss < n_seq; ss++){
                        type[ss] = pair[S[ss][i]][S[ss][j]];
                        if (type[ss]==0) type[ss] = 7;
                        en += E_ExtLoop(type[ss], -1, -1, P);
                      }
                      if (fij == en) traced=j;
                    }

                    if(with_gquad){
                      if(fij == f5[i-1] + ggg[indx[j]+i]){
                        /* found the decomposition */
                        traced = j; jj = i - 1; gq = 1;
                        break;
                      }
                    }

                    if (traced) break;
                  }
                  break;
        default:  /* j or j-1 is paired. Find pairing partner */
                  for (i=j-TURN-1,traced=0; i>=1; i--) {
                    int cc, en;
                    jj = i-1;
                    if (c[indx[j]+i]<INF) {
                      en = c[indx[j]+i] + f5[i-1];
                      for(ss = 0; ss < n_seq; ss++){
                        type[ss] = pair[S[ss][i]][S[ss][j]];
                        if (type[ss]==0) type[ss] = 7;
                        en += E_ExtLoop(type[ss], (i>1) ? S5[ss][i]: -1, (j < length) ? S3[ss][j] : -1, P);
                      }
                      if (fij == en) traced=j;
                    }

                    if(with_gquad){
                      if(fij == f5[i-1] + ggg[indx[j]+i]){
                        /* found the decomposition */
                        traced = j; jj = i - 1; gq = 1;
                        break;
                      }
                    }

                    if (traced) break;
                  }
                  break;
      }

      if (!traced) nrerror("backtrack failed in f5");
      /* push back the remaining f5 portion */
      sector[++s].i = 1;
      sector[s].j   = jj;
      sector[s].ml  = ml;

      /* trace back the base pair found */
      j=traced;

      if(with_gquad && gq){
        /* goto backtrace of gquadruplex */
        goto repeat_gquad;
      }

      base_pair2[++b].i = i;
      base_pair2[b].j   = j;
      cov_en += pscore[indx[j]+i];
      goto repeat1;
    }
    else { /* trace back in fML array */
      if (fML[indx[j]+i+1]+n_seq*P->MLbase == fij) { /* 5' end is unpaired */
        sector[++s].i = i+1;
        sector[s].j   = j;
        sector[s].ml  = ml;
        continue;
      }

      if(with_gquad){
        if(fij == ggg[indx[j]+i] + n_seq * E_MLstem(0, -1, -1, P)){
          /* go to backtracing of quadruplex */
          goto repeat_gquad;
        }
      }

      cij = c[indx[j]+i];
      if(dangles){
        for(ss = 0; ss < n_seq; ss++){
          tt = pair[S[ss][i]][S[ss][j]];
          if(tt==0) tt=7;
          cij += E_MLstem(tt, S5[ss][i], S3[ss][j], P);
        }
      }
      else{
        for(ss = 0; ss < n_seq; ss++){
          tt = pair[S[ss][i]][S[ss][j]];
          if(tt==0) tt=7;
          cij += E_MLstem(tt, -1, -1, P);
        }
      }

      if (fij==cij){
        /* found a pair */
        base_pair2[++b].i = i;
        base_pair2[b].j   = j;
        cov_en += pscore[indx[j]+i];
        goto repeat1;
      }

      for (k = i+1+TURN; k <= j-2-TURN; k++)
        if (fij == (fML[indx[k]+i]+fML[indx[j]+k+1]))
          break;

      sector[++s].i = i;
      sector[s].j   = k;
      sector[s].ml  = ml;
      sector[++s].i = k+1;
      sector[s].j   = j;
      sector[s].ml  = ml;

      if (k>j-2-TURN) nrerror("backtrack failed in fML");
      continue;
    }

  repeat1:

    /*----- begin of "repeat:" -----*/
    if (canonical)  cij = c[indx[j]+i];

    for (ss=0; ss<n_seq; ss++) {
      type[ss] = pair[S[ss][i]][S[ss][j]];
      if (type[ss]==0) type[ss] = 7;
    }

    if (noLonelyPairs)
      if (cij == c[indx[j]+i]) {
        /* (i.j) closes canonical structures, thus
           (i+1.j-1) must be a pair                */
        for (ss=0; ss<n_seq; ss++) {
          type_2 = pair[S[ss][j-1]][S[ss][i+1]];  /* j,i not i,j */
          if (type_2==0) type_2 = 7;
          cij -= P->stack[type[ss]][type_2];
        }
        cij += pscore[indx[j]+i];
        base_pair2[++b].i = i+1;
        base_pair2[b].j   = j-1;
        cov_en += pscore[indx[j-1]+i+1];
        i++; j--;
        canonical=0;
        goto repeat1;
      }
    canonical = 1;
    cij += pscore[indx[j]+i];

    {int cc=0;
    for (ss=0; ss<n_seq; ss++) {
        if ((a2s[ss][j-1]-a2s[ss][i])<3) cc+=600;
        else cc += E_Hairpin(a2s[ss][j-1]-a2s[ss][i], type[ss], S3[ss][i], S5[ss][j], Ss[ss]+a2s[ss][i-1], P);
      }
    if (cij == cc) /* found hairpin */
      continue;
    }
    for (p = i+1; p <= MIN2(j-2-TURN,i+MAXLOOP+1); p++) {
      minq = j-i+p-MAXLOOP-2;
      if (minq<p+1+TURN) minq = p+1+TURN;
      for (q = j-1; q >= minq; q--) {

        if (c[indx[q]+p]>=INF) continue;

        for (ss=energy=0; ss<n_seq; ss++) {
          type_2 = pair[S[ss][q]][S[ss][p]];  /* q,p not p,q */
          if (type_2==0) type_2 = 7;
          energy += E_IntLoop(a2s[ss][p-1]-a2s[ss][i],a2s[ss][j-1]-a2s[ss][q],
                               type[ss], type_2,
                               S3[ss][i], S5[ss][j],
                               S5[ss][p], S3[ss][q], P);

        }
        traced = (cij == energy+c[indx[q]+p]);
        if (traced) {
          base_pair2[++b].i = p;
          base_pair2[b].j   = q;
          cov_en += pscore[indx[q]+p];
          i = p, j = q;
          goto repeat1;
        }
      }
    }

    /* end of repeat: --------------------------------------------------*/

    /* (i.j) must close a multi-loop */

    i1 = i+1;
    j1 = j-1;

    if(with_gquad){
      /*
        The case that is handled here actually resembles something like
        an interior loop where the enclosing base pair is of regular
        kind and the enclosed pair is not a canonical one but a g-quadruplex
        that should then be decomposed further...
      */
      mm = 0;
      for(s=0;s<n_seq;s++){
        tt = type[s];
        if(tt == 0) tt = 7;
        if(dangles == 2)
          mm += P->mismatchI[tt][S3[s][i]][S5[s][j]];
        if(tt > 2)
          mm += P->TerminalAU;
      }

      for(p = i + 2;
        p < j - VRNA_GQUAD_MIN_BOX_SIZE;
        p++){
        if(S_cons[p] != 3) continue;
        l1    = p - i - 1;
        if(l1>MAXLOOP) break;
        minq  = j - i + p - MAXLOOP - 2;
        c0    = p + VRNA_GQUAD_MIN_BOX_SIZE - 1;
        minq  = MAX2(c0, minq);
        c0    = j - 1;
        maxq  = p + VRNA_GQUAD_MAX_BOX_SIZE + 1;
        maxq  = MIN2(c0, maxq);
        for(q = minq; q < maxq; q++){
          if(S_cons[q] != 3) continue;
          c0  = mm + ggg[indx[q] + p] + n_seq * P->internal_loop[l1 + j - q - 1];
          if(cij == c0){
            i=p;j=q;
            goto repeat_gquad;
          }
        }
      }
      p = i1;
      if(S_cons[p] == 3){
        if(p < j - VRNA_GQUAD_MIN_BOX_SIZE){
          minq  = j - i + p - MAXLOOP - 2;
          c0    = p + VRNA_GQUAD_MIN_BOX_SIZE - 1;
          minq  = MAX2(c0, minq);
          c0    = j - 3;
          maxq  = p + VRNA_GQUAD_MAX_BOX_SIZE + 1;
          maxq  = MIN2(c0, maxq);
          for(q = minq; q < maxq; q++){
            if(S_cons[q] != 3) continue;
            if(cij == mm + ggg[indx[q] + p] + n_seq * P->internal_loop[j - q - 1]){
              i = p; j=q;
              goto repeat_gquad;
            }
          }
        }
      }
      q = j1;
      if(S_cons[q] == 3)
        for(p = i1 + 3; p < j - VRNA_GQUAD_MIN_BOX_SIZE; p++){
          l1    = p - i - 1;
          if(l1>MAXLOOP) break;
          if(S_cons[p] != 3) continue;
          if(cij == mm + ggg[indx[q] + p] + n_seq * P->internal_loop[l1]){
            i = p; j = q;
            goto repeat_gquad;
          }
        }
    }

    mm = n_seq*P->MLclosing;
    if(dangles){
      for(ss = 0; ss < n_seq; ss++){
        tt = rtype[type[ss]];
        mm += E_MLstem(tt, S5[ss][j], S3[ss][i], P);
      }
    }
    else{
      for(ss = 0; ss < n_seq; ss++){
        tt = rtype[type[ss]];
        mm += E_MLstem(tt, -1, -1, P);
      }
    }
    sector[s+1].ml  = sector[s+2].ml = 1;

    for (k = i1+TURN+1; k < j1-TURN-1; k++){
      if(cij == fML[indx[k]+i1] + fML[indx[j1]+k+1] + mm) break;
    }

    if (k<=j-3-TURN) { /* found the decomposition */
      sector[++s].i = i1;
      sector[s].j   = k;
      sector[++s].i = k+1;
      sector[s].j   = j1;
    } else {
        nrerror("backtracking failed in repeat");
    }

    continue; /* this is a workarround to not accidentally proceed in the following block */

  repeat_gquad:
    /*
      now we do some fancy stuff to backtrace the stacksize and linker lengths
      of the g-quadruplex that should reside within position i,j
    */
    {
      int cnt1, cnt2, cnt3, cnt4, l[3], L, size;
      size = j-i+1;

      for(L=0; L < VRNA_GQUAD_MIN_STACK_SIZE;L++){
        if(S_cons[i+L] != 3) break;
        if(S_cons[j-L] != 3) break;
      }

      if(L == VRNA_GQUAD_MIN_STACK_SIZE){
        /* continue only if minimum stack size starting from i is possible */
        for(; L<=VRNA_GQUAD_MAX_STACK_SIZE;L++){
          if(S_cons[i+L-1] != 3) break; /* break if no more consecutive G's 5' */
          if(S_cons[j-L+1] != 3) break; /* break if no more consecutive G'1 3' */
          for(    l[0] = VRNA_GQUAD_MIN_LINKER_LENGTH;
                  (l[0] <= VRNA_GQUAD_MAX_LINKER_LENGTH)
              &&  (size - 4*L - 2*VRNA_GQUAD_MIN_LINKER_LENGTH - l[0] >= 0);
              l[0]++){
            /* check whether we find the second stretch of consecutive G's */
            for(cnt1 = 0; (cnt1 < L) && (S_cons[i+L+l[0]+cnt1] == 3); cnt1++);
            if(cnt1 < L) continue;
            for(    l[1] = VRNA_GQUAD_MIN_LINKER_LENGTH;
                    (l[1] <= VRNA_GQUAD_MAX_LINKER_LENGTH)
                &&  (size - 4*L - VRNA_GQUAD_MIN_LINKER_LENGTH - l[0] - l[1] >= 0);
                l[1]++){
              /* check whether we find the third stretch of consectutive G's */
              for(cnt1 = 0; (cnt1 < L) && (S_cons[i+2*L+l[0]+l[1]+cnt1] == 3); cnt1++);
              if(cnt1 < L) continue;

              /*
                the length of the third linker now depends on position j as well
                as the other linker lengths... so we do not have to loop too much
              */
              l[2] = size - 4*L - l[0] - l[1];
              if(l[2] < VRNA_GQUAD_MIN_LINKER_LENGTH) break;
              if(l[2] > VRNA_GQUAD_MAX_LINKER_LENGTH) continue;
              /* check for contribution */
              if(ggg[indx[j]+i] == E_gquad_ali(i, L, l, (const short **)S, n_seq, P)){
                int a;
                /* fill the G's of the quadruplex into base_pair2 */
                for(a=0;a<L;a++){
                  base_pair2[++b].i = i+a;
                  base_pair2[b].j   = i+a;
                  base_pair2[++b].i = i+L+l[0]+a;
                  base_pair2[b].j   = i+L+l[0]+a;
                  base_pair2[++b].i = i+L+l[0]+L+l[1]+a;
                  base_pair2[b].j   = i+L+l[0]+L+l[1]+a;
                  base_pair2[++b].i = i+L+l[0]+L+l[1]+L+l[2]+a;
                  base_pair2[b].j   = i+L+l[0]+L+l[1]+L+l[2]+a;
                }
                goto repeat_gquad_exit;
              }
            }
          }
        }
      }
      nrerror("backtracking failed in repeat_gquad");
    }
  repeat_gquad_exit:
    asm("nop");

  }

  /* fprintf(stderr, "covariance energy %6.2f\n", cov_en/100.);  */

  base_pair2[0].i = b;    /* save the total number of base pairs */
  free(type);
}


PUBLIC void encode_ali_sequence(const char *sequence, short *S, short *s5, short *s3, char *ss, unsigned short *as, int circular){
  unsigned int i,l;
  unsigned short p;
  l     = strlen(sequence);
  S[0]  = (short) l;
  s5[0] = s5[1] = 0;

  /* make numerical encoding of sequence */
  for(i=1; i<=l; i++){
    short ctemp;
    ctemp=(short) encode_char(toupper(sequence[i-1]));
    S[i]= ctemp ;
  }

  if (oldAliEn){
    /* use alignment sequences in all energy evaluations */
    ss[0]=sequence[0];
    for(i=1; i<l; i++){
      s5[i] = S[i-1];
      s3[i] = S[i+1];
      ss[i] = sequence[i];
      as[i] = i;
    }
    ss[l]   = sequence[l];
    as[l]   = l;
    s5[l]   = S[l-1];
    s3[l]   = 0;
    S[l+1]  = S[1];
    s5[1]   = 0;
    if (circular) {
      s5[1]   = S[l];
      s3[l]   = S[1];
      ss[l+1] = S[1];
    }
  }
  else{
    if(circular){
      for(i=l; i>0; i--){
        char c5;
        c5 = sequence[i-1];
        if ((c5=='-')||(c5=='_')||(c5=='~')||(c5=='.')) continue;
        s5[1] = S[i];
        break;
      }
      for (i=1; i<=l; i++) {
        char c3;
        c3 = sequence[i-1];
        if ((c3=='-')||(c3=='_')||(c3=='~')||(c3=='.')) continue;
        s3[l] = S[i];
        break;
      }
    }
    else  s5[1]=s3[l]=0;

    for(i=1,p=0; i<=l; i++){
      char c5;
      c5 = sequence[i-1];
      if ((c5=='-')||(c5=='_')||(c5=='~')||(c5=='.'))
        s5[i+1]=s5[i];
      else { /* no gap */
        ss[p++]=sequence[i-1]; /*start at 0!!*/
        s5[i+1]=S[i];
      }
      as[i]=p;
    }
    for (i=l; i>=1; i--) {
      char c3;
      c3 = sequence[i-1];
      if ((c3=='-')||(c3=='_')||(c3=='~')||(c3=='.'))
        s3[i-1]=s3[i];
      else
        s3[i-1]=S[i];
    }
  }
}

PRIVATE void make_pscores(const short *const* S, const char **AS,
                          int n_seq, const char *structure) {
  /* calculate co-variance bonus for each pair depending on  */
  /* compensatory/consistent mutations and incompatible seqs */
  /* should be 0 for conserved pairs, >0 for good pairs      */
#define NONE -10000 /* score for forbidden pairs */
  int n,i,j,k,l,s;

  int olddm[7][7]={{0,0,0,0,0,0,0}, /* hamming distance between pairs */
                  {0,0,2,2,1,2,2} /* CG */,
                  {0,2,0,1,2,2,2} /* GC */,
                  {0,2,1,0,2,1,2} /* GU */,
                  {0,1,2,2,0,2,1} /* UG */,
                  {0,2,2,1,2,0,2} /* AU */,
                  {0,2,2,2,1,2,0} /* UA */};

  float **dm;
  n=S[0][0];  /* length of seqs */
  if (ribo) {
    if (RibosumFile !=NULL) dm=readribosum(RibosumFile);
    else dm=get_ribosum(AS,n_seq,n);
  }
  else { /*use usual matrix*/
    dm=(float **)space(7*sizeof(float*));
    for (i=0; i<7;i++) {
      dm[i]=(float *)space(7*sizeof(float));
      for (j=0; j<7; j++)
        dm[i][j] = (float) olddm[i][j];
    }
  }

  for (i=1; i<n; i++) {
    for (j=i+1; (j<i+TURN+1) && (j<=n); j++)
      pscore[indx[j]+i] = NONE;
    for (j=i+TURN+1; j<=n; j++) {
      int pfreq[8]={0,0,0,0,0,0,0,0};
      double score;
      for (s=0; s<n_seq; s++) {
        int type;
        if (S[s][i]==0 && S[s][j]==0) type = 7; /* gap-gap  */
        else {
          if ((AS[s][i] == '~')||(AS[s][j] == '~')) type = 7;
          else type = pair[S[s][i]][S[s][j]];
        }
        pfreq[type]++;
      }
      if (pfreq[0]*2+pfreq[7]>n_seq) { pscore[indx[j]+i] = NONE; continue;}
      for (k=1,score=0; k<=6; k++) /* ignore pairtype 7 (gap-gap) */
        for (l=k; l<=6; l++)
          score += pfreq[k]*pfreq[l]*dm[k][l];
      /* counter examples score -1, gap-gap scores -0.25   */
      pscore[indx[j]+i] = cv_fact *
        ((UNIT*score)/n_seq - nc_fact*UNIT*(pfreq[0] + pfreq[7]*0.25));
    }
  }

  if (noLonelyPairs) /* remove unwanted pairs */
    for (k=1; k<n-TURN-1; k++)
      for (l=1; l<=2; l++) {
        int type,ntype=0,otype=0;
        i=k; j = i+TURN+l;
        type = pscore[indx[j]+i];
        while ((i>=1)&&(j<=n)) {
          if ((i>1)&&(j<n)) ntype = pscore[indx[j+1]+i-1];
          if ((otype<cv_fact*MINPSCORE)&&(ntype<cv_fact*MINPSCORE))  /* too many counterexamples */
            pscore[indx[j]+i] = NONE; /* i.j can only form isolated pairs */
          otype =  type;
          type  = ntype;
          i--; j++;
        }
      }


  if (fold_constrained&&(structure!=NULL)) {
    int psij, hx, hx2, *stack, *stack2;
    stack = (int *) space(sizeof(int)*(n+1));
    stack2 = (int *) space(sizeof(int)*(n+1));

    for(hx=hx2=0, j=1; j<=n; j++) {
      switch (structure[j-1]) {
      case 'x': /* can't pair */
        for (l=1; l<j-TURN; l++) pscore[indx[j]+l] = NONE;
        for (l=j+TURN+1; l<=n; l++) pscore[indx[l]+j] = NONE;
        break;
      case '(':
        stack[hx++]=j;
        /* fallthrough */
      case '[':
        stack2[hx2++]=j;
        /* fallthrough */
      case '<': /* pairs upstream */
        for (l=1; l<j-TURN; l++) pscore[indx[j]+l] = NONE;
        break;
      case ']':
        if (hx2<=0) {
          fprintf(stderr, "%s\n", structure);
          nrerror("unbalanced brackets in constraints");
        }
        i = stack2[--hx2];
        pscore[indx[j]+i]=NONE;
        break;
      case ')':
        if (hx<=0) {
          fprintf(stderr, "%s\n", structure);
          nrerror("unbalanced brackets in constraints");
        }
        i = stack[--hx];
        psij = pscore[indx[j]+i]; /* store for later */
        for (k=j; k<=n; k++)
          for (l=i; l<=j; l++)
            pscore[indx[k]+l] = NONE;
        for (l=i; l<=j; l++)
          for (k=1; k<=i; k++)
            pscore[indx[l]+k] = NONE;
        for (k=i+1; k<j; k++)
          pscore[indx[k]+i] = pscore[indx[j]+k] = NONE;
        pscore[indx[j]+i] = (psij>0) ? psij : 0;
        /* fallthrough */
      case '>': /* pairs downstream */
        for (l=j+TURN+1; l<=n; l++) pscore[indx[l]+j] = NONE;
        break;
      }
    }
    if (hx!=0) {
      fprintf(stderr, "%s\n", structure);
      nrerror("unbalanced brackets in constraint string");
    }
    free(stack); free(stack2);
  }
  /*free dm */
  for (i=0; i<7;i++) {
    free(dm[i]);
  }
  free(dm);
}

/*--------New scoring part-----------------------------------*/

PUBLIC int get_mpi(char *Alseq[], int n_seq, int length, int *mini) {
  int i, j,k;
  float ident=0;
  int pairnum=0;
  int sumident=0;
  float minimum=1.;
  for(j=0; j<n_seq-1; j++)
    for(k=j+1; k<n_seq; k++) {
      ident=0;
      for (i=1; i<=length; i++){
        if (Alseq[k][i]==Alseq[j][i]) ident++;
        pairnum++;
      }
      if ((ident/length)<minimum) minimum=ident/(float)length;
      sumident+=ident;
    }
  mini[0]=(int)(minimum*100.);
  if (pairnum>0)   return (int) (sumident*100/pairnum);
  else return 0;

}

PUBLIC float **readribosum(char *name){

  float **dm;
  char *line;
  FILE *fp;
  int i=0;
  int who=0;
  float a,b,c,d,e,f;
  int translator[7]={0,5,1,2,3,6,4};

  fp=fopen(name,"r");
  dm=(float **)space(7*sizeof(float*));
  for (i=0; i<7;i++) {
    dm[i]=(float *)space(7*sizeof(float));
  }
  while(1) { /*bisma hoit fertisch san*/
    line=get_line(fp);
    if (*line=='#') continue;
    i=0;
    i=sscanf(line,"%f %f %f %f %f %f",&a,&b,&c,&d,&e,&f);
    if (i==0) break;
    dm[translator[++who]][translator[1]]=a;
    dm[translator[who]][translator[2]]=b;
    dm[translator[who]][translator[3]]=c;
    dm[translator[who]][translator[4]]=d;
    dm[translator[who]][translator[5]]=e;
    dm[translator[who]][translator[6]]=f;
    free(line);
    if (who==6) break;
  }
  fclose(fp);
  return dm;
}


PRIVATE void en_corr_of_loop_gquad(int i,
                                  int j,
                                  const char **sequences,
                                  const char *structure,
                                  short *pt,
                                  int *loop_idx,
                                  int n_seq,
                                  int en[2]){

  int pos, energy, en_covar, p, q, r, s, u, type, type2, gq_en[2];
  int L, l[3];

  energy = en_covar = 0;
  q = i;
  while((pos = parse_gquad(structure + q-1, &L, l)) > 0){
    q += pos-1;
    p = q - 4*L - l[0] - l[1] - l[2] + 1;
    if(q > j) break;
    /* we've found the first g-quadruplex at position [p,q] */
    E_gquad_ali_en(p, L, l, (const short **)S, n_seq, gq_en, P);
    energy    += gq_en[0];
    en_covar  += gq_en[1];
    /* check if it's enclosed in a base pair */
    if(loop_idx[p] == 0){ q++; continue; /* g-quad in exterior loop */}
    else{
      energy += E_MLstem(0, -1, -1, P) * n_seq;
      /*  find its enclosing pair */
      int num_elem, num_g, elem_i, elem_j, up_mis;
      num_elem  = 0;
      num_g     = 1;
      r         = p - 1;
      up_mis    = q - p + 1;

      /* seek for first pairing base located 5' of the g-quad */
      for(r = p - 1; !pt[r] && (r >= i); r--);
      if(r < i) nrerror("this should not happen");

      if(r < pt[r]){ /* found the enclosing pair */
        s = pt[r];
      } else {
        num_elem++;
        elem_i = pt[r];
        elem_j = r;
        r = pt[r]-1 ;
        /* seek for next pairing base 5' of r */
        for(; !pt[r] && (r >= i); r--);
        if(r < i) nrerror("so nich");
        if(r < pt[r]){ /* found the enclosing pair */
          s = pt[r];
        } else {
          /* hop over stems and unpaired nucleotides */
          while((r > pt[r]) && (r >= i)){
            if(pt[r]){ r = pt[r]; num_elem++;}
            r--;
          }
          if(r < i) nrerror("so nich");
          s = pt[r]; /* found the enclosing pair */
        }
      }
      /* now we have the enclosing pair (r,s) */

      u = q+1;
      /* we know everything about the 5' part of this loop so check the 3' part */
      while(u<s){
        if(structure[u-1] == '.') u++;
        else if (structure[u-1] == '+'){ /* found another gquad */
          pos = parse_gquad(structure + u - 1, &L, l);
          if(pos > 0){
            E_gquad_ali_en(u, L, l, (const short **)S, n_seq, gq_en, P);
            energy += gq_en[0] + E_MLstem(0, -1, -1, P) * n_seq;
            en_covar += gq_en[1];
            up_mis += pos;
            u += pos;
            num_g++;
          }
        } else { /* we must have found a stem */
          if(!(u < pt[u])) nrerror("wtf!");
          num_elem++;
          elem_i = u;
          elem_j = pt[u];
          en_corr_of_loop_gquad(u,
                                pt[u],
                                sequences,
                                structure,
                                pt,
                                loop_idx,
                                n_seq,
                                gq_en);
          energy    += gq_en[0];
          en_covar  += gq_en[1];
          u = pt[u] + 1;
        }
      }
      if(u!=s) nrerror("what the ...");
      else{ /* we are done since we've found no other 3' structure element */
        switch(num_elem){
          /* g-quad was misinterpreted as hairpin closed by (r,s) */
          case 0:   /*if(num_g == 1)
                      if((p-r-1 == 0) || (s-q-1 == 0))
                        nrerror("too few unpaired bases");
                    */
                    {
                      int ee = 0;
                      int cnt;
                      for(cnt=0;cnt<n_seq;cnt++){
                        type = pair[S[cnt][r]][S[cnt][s]];
                        if(type == 0) type = 7;
                        if ((a2s[cnt][s-1]-a2s[cnt][r])<3) ee+=600;
                        else ee += E_Hairpin( a2s[cnt][s-1] - a2s[cnt][r],
                                              type,
                                              S3[cnt][r],
                                              S5[cnt][s],
                                              Ss[cnt] + a2s[cnt][r-1],
                                              P);
                      }
                      energy -= ee;
                      ee = 0;
                      for(cnt=0;cnt < n_seq; cnt++){
                        type = pair[S[cnt][r]][S[cnt][s]];
                        if(type == 0) type = 7;
                        if(dangles == 2)
                          ee += P->mismatchI[type][S3[cnt][r]][S5[cnt][s]];
                        if(type > 2)
                          ee += P->TerminalAU;
                      }
                      energy += ee;
                    }
                    energy += n_seq * P->internal_loop[s-r-1-up_mis];
                    break;
          /* g-quad was misinterpreted as interior loop closed by (r,s) with enclosed pair (elem_i, elem_j) */
          case 1:   {
                      int ee = 0;
                      int cnt;
                      for(cnt = 0; cnt<n_seq;cnt++){
                        type = pair[S[cnt][r]][S[cnt][s]];
                        if(type == 0) type = 7;
                        type2 = pair[S[cnt][elem_j]][S[cnt][elem_i]];
                        if(type2 == 0) type2 = 7;
                        ee += E_IntLoop(a2s[cnt][elem_i-1] - a2s[cnt][r],
                                        a2s[cnt][s-1] - a2s[cnt][elem_j],
                                        type,
                                        type2,
                                        S3[cnt][r],
                                        S5[cnt][s],
                                        S5[cnt][elem_i],
                                        S3[cnt][elem_j],
                                        P);
                      }
                      energy -= ee;
                      ee = 0;
                      for(cnt = 0; cnt < n_seq; cnt++){
                        type = pair[S[cnt][s]][S[cnt][r]];
                        if(type == 0) type = 7;
                        ee += E_MLstem(type, S5[cnt][s], S3[cnt][r], P);
                        type = pair[S[cnt][elem_i]][S[cnt][elem_j]];
                        if(type == 0) type = 7;
                        ee += E_MLstem(type, S5[cnt][elem_i], S3[cnt][elem_j], P);
                      }
                      energy += ee;
                    }
                    energy += (P->MLclosing + (elem_i-r-1+s-elem_j-1-up_mis) * P->MLbase) * n_seq;
                    break;
          /* gquad was misinterpreted as unpaired nucleotides in a multiloop */
          default:  energy -= (up_mis) * P->MLbase * n_seq;
                    break;
        }
      }
      q = s+1;
    }
  }
  en[0] = energy;
  en[1] = en_covar;
}

PUBLIC float energy_of_ali_gquad_structure( const char **sequences,
                                            const char *structure,
                                            int n_seq,
                                            float *energy){

  int new=0;
  unsigned int s, n;
  short *pt;

  short           **tempS;
  short           **tempS5;     /*S5[s][i] holds next base 5' of i in sequence s*/
  short           **tempS3;     /*Sl[s][i] holds next base 3' of i in sequence s*/
  char            **tempSs;
  unsigned short  **tempa2s;

  int *tempindx, *loop_idx;
  int *temppscore;

  int en_struct[2], gge[2];

  if(sequences[0] != NULL){
    n = (unsigned int) strlen(sequences[0]);
    update_alifold_params();

    /*save old memory*/
    tempS = S; tempS3 = S3; tempS5 = S5; tempSs = Ss; tempa2s = a2s;
    tempindx = indx; temppscore = pscore;

    alloc_sequence_arrays(sequences, &S, &S5, &S3, &a2s, &Ss, 0);
    pscore  = (int *) space(sizeof(int)*((n+1)*(n+2)/2));
    indx    = get_indx(n);
    make_pscores((const short *const*)S, sequences, n_seq, NULL);
    new     = 1;

    pt = make_pair_table(structure);
    energy_of_alistruct_pt(sequences,pt, n_seq, &(en_struct[0]));
    loop_idx    = make_loop_index_pt(pt);
    en_corr_of_loop_gquad(1, n, sequences, structure, pt, loop_idx, n_seq, gge);
    en_struct[0] += gge[0];
    en_struct[1] += gge[1];

    free(loop_idx);
    free(pt);
    energy[0] = (float)en_struct[0]/(float)(100*n_seq);
    energy[1] = (float)en_struct[1]/(float)(100*n_seq);

    free(pscore);
    free(indx);
    free_sequence_arrays(n_seq, &S, &S5, &S3, &a2s, &Ss);

    /* restore old memory */
    S = tempS; S3 = tempS3; S5 = tempS5; Ss = tempSs; a2s = tempa2s;
    indx = tempindx; pscore = temppscore;
  }
  else nrerror("energy_of_alistruct(): no sequences in alignment!");

  return energy[0];

}

PUBLIC  float energy_of_alistruct(const char **sequences, const char *structure, int n_seq, float *energy){
  int new=0;
  unsigned int s, n;
  short *pt;

  short           **tempS;
  short           **tempS5;     /*S5[s][i] holds next base 5' of i in sequence s*/
  short           **tempS3;     /*Sl[s][i] holds next base 3' of i in sequence s*/
  char            **tempSs;
  unsigned short  **tempa2s;

  int *tempindx;
  int *temppscore;

  int en_struct[2];

  if(sequences[0] != NULL){
    n = (unsigned int) strlen(sequences[0]);
    update_alifold_params();

    /*save old memory*/
    tempS = S; tempS3 = S3; tempS5 = S5; tempSs = Ss; tempa2s = a2s;
    tempindx = indx; temppscore = pscore;

    alloc_sequence_arrays(sequences, &S, &S5, &S3, &a2s, &Ss, 0);
    pscore  = (int *) space(sizeof(int)*((n+1)*(n+2)/2));
    indx    = get_indx(n);
    make_pscores((const short *const*)S, sequences, n_seq, NULL);
    new     = 1;

    pt = make_pair_table(structure);
    energy_of_alistruct_pt(sequences,pt, n_seq, &(en_struct[0]));

    free(pt);
    energy[0] = (float)en_struct[0]/(float)(100*n_seq);
    energy[1] = (float)en_struct[1]/(float)(100*n_seq);

    free(pscore);
    free(indx);
    free_sequence_arrays(n_seq, &S, &S5, &S3, &a2s, &Ss);

    /* restore old memory */
    S = tempS; S3 = tempS3; S5 = tempS5; Ss = tempSs; a2s = tempa2s;
    indx = tempindx; pscore = temppscore;
  }
  else nrerror("energy_of_alistruct(): no sequences in alignment!");

  return energy[0];
}

PRIVATE void energy_of_alistruct_pt(const char **sequences, short *pt, int n_seq, int *energy){
  int i, length;

  length = S[0][0];
  energy[0] =  backtrack_type=='M' ? ML_Energy_pt(0, n_seq, pt) : EL_Energy_pt(0, n_seq, pt);
  energy[1] = 0;
  for (i=1; i<=length; i++) {
    if (pt[i]==0) continue;
    stack_energy_pt(i, sequences, pt, n_seq, energy);
    i=pt[i];
  }
}

PRIVATE void stack_energy_pt(int i, const char **sequences, short *pt, int n_seq, int *energy)
{
  /* calculate energy of substructure enclosed by (i,j) */
  int ee= 0;
  int j, p, q, s;
  int *type = (int *) space(n_seq*sizeof(int));

  j = pt[i];
  for (s=0; s<n_seq; s++) {
    type[s] = pair[S[s][i]][S[s][j]];
    if (type[s]==0) {
    type[s]=7;
    }
  }
  p=i; q=j;
  while (p<q) { /* process all stacks and interior loops */
    int type_2;
    while (pt[++p]==0);
    while (pt[--q]==0);
    if ((pt[q]!=(short)p)||(p>q)) break;
    ee=0;
    for (s=0; s<n_seq; s++) {
      type_2 = pair[S[s][q]][S[s][p]];
      if (type_2==0) {
        type_2=7;
      }
      ee += E_IntLoop(a2s[s][p-1]-a2s[s][i], a2s[s][j-1]-a2s[s][q], type[s], type_2, S3[s][i], S5[s][j], S5[s][p], S3[s][q], P);
    }
    energy[0] += ee;
    energy[1] += pscore[indx[j]+i];
    i=p; j=q;
    for (s=0; s<n_seq; s++) {
      type[s] = pair[S[s][i]][S[s][j]];
      if (type[s]==0) type[s]=7;
    }
  }  /* end while */

  /* p,q don't pair must have found hairpin or multiloop */

  if (p>q) {
    ee=0;/* hair pin */
    for (s=0; s< n_seq; s++) {
      if ((a2s[s][j-1]-a2s[s][i])<3) ee+=600;
      else ee += E_Hairpin(a2s[s][j-1]-a2s[s][i], type[s], S3[s][i], S5[s][j], Ss[s]+(a2s[s][i-1]), P);
    }
    energy[0] += ee;
    energy[1] += pscore[indx[j]+i];
    free(type);
    return;
  }
  /* (i,j) is exterior pair of multiloop */
  energy[1] += pscore[indx[j]+i];
  while (p<j) {
    /* add up the contributions of the substructures of the ML */
    stack_energy_pt(p, sequences, pt, n_seq, energy);
    p = pt[p];
    /* search for next base pair in multiloop */
    while (pt[++p]==0);
  }
  energy[0] += ML_Energy_pt(i, n_seq, pt);
  free(type);
}



PRIVATE int ML_Energy_pt(int i, int n_seq, short *pt){
  /* i is the 5'-base of the closing pair */

  int   energy, tt, i1, j, p, q, u, s;
  short d5, d3;

  j = pt[i];
  i1  = i;
  p   = i+1;
  u   = 0;
  energy = 0;

  do{ /* walk around the multi-loop */
    /* hop over unpaired positions */
    while (p < j && pt[p]==0) p++;
    if(p >= j) break;
    /* get position of pairing partner */
    q  = pt[p];
    /* memorize number of unpaired positions? no, we just approximate here */
    u += p-i1-1;

    for (s=0; s< n_seq; s++) {
      /* get type of base pair P->q */
      tt = pair[S[s][p]][S[s][q]];
      if (tt==0) tt=7;
      d5 = dangles && (a2s[s][p]>1) && (tt!=0) ? S5[s][p] : -1;
      d3 = dangles && (a2s[s][q]<a2s[s][S[0][0]]) ? S3[s][q] : -1;
      energy += E_MLstem(tt, d5, d3, P);
    }
    i1  = q;
    p   = q + 1;
  }while(1);

  if(i > 0){
    energy  += P->MLclosing * n_seq;
    if(dangles){
      for (s=0; s<n_seq; s++){
        tt = pair[S[s][j]][S[s][i]];
        if (tt==0) tt=7;
        energy += E_MLstem(tt, S5[s][j], S3[s][i], P);
      }
    }
    else{
      for (s=0; s<n_seq; s++){
        tt = pair[S[s][j]][S[s][i]];
        if (tt==0) tt=7;
        energy += E_MLstem(tt, -1, -1, P);
      }
    }
  }
  u += j - i1 - 1;
  energy += u * P->MLbase * n_seq;
  return energy;
}

PRIVATE int EL_Energy_pt(int i, int n_seq, short *pt){
  int   energy, tt, i1, j, p, q, s;
  short d5, d3;

  j = pt[0];

  p   = i+1;
  energy = 0;

  do{ /* walk along the backbone */
    /* hop over unpaired positions */
    while (p < j && pt[p]==0) p++;
    if(p == j) break; /* no more stems */
    /* get position of pairing partner */
    q  = pt[p];
    for (s=0; s< n_seq; s++) {
      /* get type of base pair P->q */
      tt = pair[S[s][p]][S[s][q]];
      if (tt==0) tt=7;
      d5 = dangles && (a2s[s][p]>1) && (tt!=0) ? S5[s][p] : -1;
      d3 = dangles && (a2s[s][q]<a2s[s][S[0][0]]) ? S3[s][q] : -1;
      energy += E_ExtLoop(tt, d5, d3, P);
    }
    p   = q + 1;
  }while(p < j);

  return energy;
}