/*
                  minimum free energy
                  RNA secondary structure prediction
                  with maximum distance base pairs

                  c Ivo Hofacker, Peter Stadler

                  Vienna RNA package
*/

#ifdef HAVE_CONFIG_H
#include "config.h"
#endif

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <ctype.h>
#include <string.h>
#include <limits.h>
#include "ViennaRNA/utils.h"
#include "ViennaRNA/energy_par.h"
#include "ViennaRNA/fold_vars.h"
#include "ViennaRNA/params.h"
#include "ViennaRNA/loop_energies.h"
#include "ViennaRNA/gquad.h"
#include "ViennaRNA/Lfold.h"

#ifdef USE_SVM
#include "svm.h"
#include "svm_utils.h"
#endif

#define MAXSECTORS                  500   /* dimension for a backtrack array */
#define INT_CLOSE_TO_UNDERFLOW(i)   ((i) <= (INT_MIN/16))
#define UNDERFLOW_CORRECTION        (INT_MIN/32)

/*
#################################
# GLOBAL VARIABLES              #
#################################
*/

/*
#################################
# PRIVATE VARIABLES             #
#################################
*/

/*
#################################
# PRIVATE FUNCTION DECLARATIONS #
#################################
*/
PRIVATE float wrap_Lfold( vrna_fold_compound_t *vc,
                          int with_zsc,
                          double min_z,
                          FILE *file);
PRIVATE void  make_ptypes(vrna_fold_compound_t *vc,
                          int i);
PRIVATE char  *backtrack( vrna_fold_compound_t *vc,
                          int start,
                          int maxdist);
PRIVATE int   fill_arrays(vrna_fold_compound_t *vc,
                          int with_zsc,
                          double min_z,
#ifdef USE_SVM
                          struct svm_model *avg_model,
                          struct svm_model *sd_model,
#endif
                          int *underflow,
                          FILE *output);

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

PUBLIC float
vrna_Lfold( const char *string,
            int window_size,
            FILE  *file){

  float               energy;
  vrna_fold_compound_t  *vc;
  vrna_md_t           md;

  vrna_md_set_default(&md);

  md.window_size = window_size;
  md.max_bp_span = window_size;

  vc  = vrna_fold_compound(string, &md, VRNA_OPTION_WINDOW);

  energy = wrap_Lfold(vc, 0, 0.0, file);

  vrna_fold_compound_free(vc);

  return energy;
}

PUBLIC float
vrna_mfe_window( vrna_fold_compound_t *vc,
            FILE *file){

  return wrap_Lfold(vc, 0, 0.0, file);
}

#ifdef USE_SVM

PUBLIC float
vrna_Lfoldz(const char *string,
            int window_size,
            double min_z,
            FILE *file){

  float               energy;
  vrna_fold_compound_t  *vc;
  vrna_md_t           md;

  vrna_md_set_default(&md);

  md.window_size = window_size;
  md.max_bp_span = window_size;

  vc  = vrna_fold_compound(string, &md, VRNA_OPTION_WINDOW);

  energy = wrap_Lfold(vc, 1, min_z, file);

  vrna_fold_compound_free(vc);

  return energy;
}

PUBLIC float
vrna_mfe_window_zscore( vrna_fold_compound_t *vc,
             double min_z,
             FILE *file){

  return wrap_Lfold(vc, 1, min_z, file);
}

#endif

/*
#####################################
# BEGIN OF STATIC HELPER FUNCTIONS  #
#####################################
*/

PRIVATE float
wrap_Lfold( vrna_fold_compound_t *vc,
            int with_zsc,
            double min_z,
            FILE *file){

  int     i, energy, underflow, n, maxdist;
  float   mfe_local;
  FILE    *out;

#ifdef USE_SVM
  struct svm_model  *avg_model = NULL;
  struct svm_model  *sd_model = NULL;
#endif

  if(!vrna_fold_compound_prepare(vc, VRNA_OPTION_MFE | VRNA_OPTION_WINDOW)){
    vrna_message_warning("vrna_mfe_window@Lfold.c: Failed to prepare vrna_fold_compound");
    return (float)(INF/100.);
  }


  n       = vc->length;
  maxdist = vc->window_size;
  out     = (file) ? file : stdout;

  for(i = n; (i >= (int)n - (int)maxdist - 4) && (i > 0); i--)
    make_ptypes(vc, i);

#ifdef USE_SVM  /*svm*/
  if(with_zsc){
    avg_model = svm_load_model_string(avg_model_string);
    sd_model  = svm_load_model_string(sd_model_string);
  }
#endif

  /* keep track of how many times we were close to an integer underflow */
  underflow = 0;

#ifdef USE_SVM
  energy = fill_arrays(vc, with_zsc, min_z, avg_model, sd_model, &underflow, out);
  if(with_zsc){
    svm_free_model_content(avg_model);
    svm_free_model_content(sd_model);
  }
#else
  energy = fill_arrays(vc, with_zsc, min_z, &underflow, out);
#endif

  mfe_local = (underflow > 0) ? ((float)underflow * (float)(UNDERFLOW_CORRECTION)) / 100. : 0.;
  mfe_local += (float)energy/100.;

  return mfe_local;
}

PRIVATE int
fill_arrays(vrna_fold_compound_t *vc,
            int zsc,
            double min_z,
#ifdef USE_SVM
            struct svm_model *avg_model,
            struct svm_model *sd_model,
#endif
            int *underflow,
            FILE *output){

  /* fill "c", "fML" and "f3" arrays and return  optimal energy */

  int   i, j, k, length, energy, maxdist;
  int   **c, **fML, *f3, **ggg;
  int   decomp, new_fML;
  int   no_close, type, type_2, tt, with_gquad, dangle_model, noLP, noGUclosure, turn;
  int   *rtype;
  int   fij;
  int   lind;

  int   *cc = NULL;        /* linear array for calculating canonical structures */
  int   *cc1 = NULL;       /*   "     "        */
  int   *Fmi = NULL;       /* holds row i of fML (avoids jumps in memory) */
  int   *DMLi = NULL;      /* DMLi[j] holds MIN(fML[i,k]+fML[k+1,j])  */
  int   *DMLi1 = NULL;     /*             MIN(fML[i+1,k]+fML[k+1,j])  */
  int   *DMLi2 = NULL;     /*             MIN(fML[i+2,k]+fML[k+1,j])  */

  short         *S, *S1;
  char          *string, **ptype, *prev;
  vrna_param_t  *P;
  vrna_md_t     *md;


  string        = vc->sequence;
  length        = vc->length;
  S             = vc->sequence_encoding2;
  S1            = vc->sequence_encoding;
  ptype         = vc->ptype_local;
  maxdist       = vc->window_size;
  P             = vc->params;
  md            = &(P->model_details);
  dangle_model  = md->dangles;
  with_gquad    = md->gquad;
  noLP          = md->noLP;
  noGUclosure   = md->noGUclosure;
  turn          = md->min_loop_size;
  rtype         = &(md->rtype[0]);

  prev          = NULL;

  c           = vc->matrices->c_local;
  fML         = vc->matrices->fML_local;
  f3          = vc->matrices->f3_local;
  ggg         = vc->matrices->ggg_local;

  cc    = (int *)   vrna_alloc(sizeof(int)   *(maxdist+5));
  cc1   = (int *)   vrna_alloc(sizeof(int)   *(maxdist+5));
  Fmi   = (int *)   vrna_alloc(sizeof(int)   *(maxdist+5));
  DMLi  = (int *)   vrna_alloc(sizeof(int)   *(maxdist+5));
  DMLi1 = (int *)   vrna_alloc(sizeof(int)   *(maxdist+5));
  DMLi2 = (int *)   vrna_alloc(sizeof(int)   *(maxdist+5));


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

  if(with_gquad){
    vrna_gquad_mx_local_update(vc, length - maxdist - 4);
    ggg = vc->matrices->ggg_local;
  }

  for (i = length-turn-1; i >= 1; i--) { /* i,j in [1..length] */
    for (j = i+turn+1; j <= length && j <= i+maxdist; j++) {
      int p, q;
      type = ptype[i][j-i];

      no_close = (((type==3)||(type==4))&&noGUclosure);

      if (type) {   /* we have a pair */
        int new_c=0, stackEnergy=INF;
        /* hairpin ----------------------------------------------*/

        new_c = (no_close) ? FORBIDDEN : E_Hairpin(j-i-1, type, S1[i+1], S1[j-1], string+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++){
          int minq = j-i+p-MAXLOOP-2;
          if (minq<p+1+turn) minq = p+1+turn;
          for (q = minq; q < j; q++) {
            type_2 = ptype[p][q-p];

            if (type_2==0) continue;
            type_2 = rtype[type_2];

            if (noGUclosure)
              if (no_close||(type_2==3)||(type_2==4))
                if ((p>i+1)||(q<j-1)) continue;  /* continue unless stack */

            energy = E_IntLoop(p-i-1, j-q-1, type, type_2, S1[i+1], S1[j-1], S1[p-1], S1[q+1],P);
            new_c = MIN2(new_c, energy + c[p][q-p]);
            if ((p==i+1)&&(j==q+1)) stackEnergy = energy; /* remember stack energy */
          } /* end q-loop */
        } /* end p-loop */

        /* multi-loop decomposition ------------------------*/
        if (!no_close) {
          decomp  = DMLi1[j-1-(i+1)];
          tt      = rtype[type];
          switch(dangle_model){
            /* no dangle_model */
            case 0:   decomp += E_MLstem(tt, -1, -1, P);
                      break;
            /* double dangle_model */
            case 2:   decomp += E_MLstem(tt, S1[j-1], S1[i+1], P);
                      break;
            /* normal dangle_model, aka dangle_model = 1 */
            default:  decomp += E_MLstem(tt, -1, -1, P);
                      decomp = MIN2(decomp, DMLi2[j-1-(i+2)] + E_MLstem(tt, -1, S1[i+1], P) + P->MLbase);
                      decomp = MIN2(decomp, DMLi2[j-2-(i+2)] + E_MLstem(tt, S1[j-1], S1[i+1], P) + 2*P->MLbase);
                      decomp = MIN2(decomp, DMLi1[j-2-(i+1)] + E_MLstem(tt, S1[j-1], -1, P) + P->MLbase);
                      break;
          }
          new_c = MIN2(new_c, decomp + P->MLclosing);
        }

        /* coaxial stacking of (i.j) with (i+1.k) or (k+1.j-1) */

        if (dangle_model==3) {
          decomp = INF;
          for (k = i+2+turn; k < j-2-turn; k++) {
            type_2 = ptype[i+1][k-i-1]; type_2 = rtype[type_2];
            if (type_2)
              decomp = MIN2(decomp, c[i+1][k-i-1]+P->stack[type][type_2]+
                            fML[k+1][j-1-k-1]);
            type_2 = ptype[k+1][j-1-k-1]; type_2 = rtype[type_2];
            if (type_2)
              decomp = MIN2(decomp, c[k+1][j-1-k-1]+P->stack[type][type_2]+
                            fML[i+1][k-i-1]);
          }
          /* no TermAU penalty if coax stack */
          decomp += 2*P->MLintern[1] + P->MLclosing;
          new_c = MIN2(new_c, decomp);
        }

        if(with_gquad){
          /* include all cases where a g-quadruplex may be enclosed by base pair (i,j) */
          if (!no_close) {
            tt = rtype[type];
            energy = E_GQuad_IntLoop_L(i, j, type, S1, ggg, maxdist, P);
            new_c = MIN2(new_c, energy);
          }
        }

        new_c = MIN2(new_c, cc1[j-1-(i+1)]+stackEnergy);
        cc[j-i] = new_c;
        if (noLP)
          c[i][j-i] = cc1[j-1-(i+1)]+stackEnergy;
        else
          c[i][j-i] = cc[j-i];

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

      else c[i][j-i] = INF;

      /* done with c[i,j], now compute fML[i,j] */
      /* free ends ? -----------------------------------------*/
      new_fML = INF;
      switch(dangle_model){
        /* no dangle_model */
        case 0:   new_fML = fML[i+1][j-i-1] + P->MLbase;
                  new_fML = MIN2(new_fML, fML[i][j-1-i] + P->MLbase);
                  new_fML = MIN2(new_fML, c[i][j-i] + E_MLstem(type, -1, -1, P));
                  break;
        /* double dangle_model */
        case 2:   new_fML = fML[i+1][j-i-1] + P->MLbase;
                  new_fML = MIN2(fML[i][j-1-i] + P->MLbase, new_fML);
                  new_fML = MIN2(new_fML,  c[i][j-i] + E_MLstem(type, (i>1) ? S1[i-1] : -1, (j<length) ? S1[j+1] : -1, P));
                  break;
        /* normal dangle_model, aka dangle_model = 1 */
        default:  /* i unpaired */
                  new_fML = fML[i+1][j-i-1] + P->MLbase;
                  /* j unpaired */
                  new_fML = MIN2(new_fML, fML[i][j-1-i] + P->MLbase);
                  /* i,j */
                  if(type) new_fML = MIN2(new_fML, c[i][j-i] + E_MLstem(type, -1, -1, P));
                  /* i+1,j */
                  tt = ptype[i+1][j-i-1];
                  if(tt) new_fML = MIN2(new_fML, c[i+1][j-i-1] + E_MLstem(tt, S1[i], -1, P) + P->MLbase);
                  /* i, j-1 */
                  tt = ptype[i][j-1-i];
                  if(tt) new_fML = MIN2(new_fML, c[i][j-1-i] + E_MLstem(tt, -1, S1[j], P) + P->MLbase);
                  /* i+1,j-1 */
                  tt = ptype[i+1][j-1-i-1];
                  if(tt) new_fML = MIN2(new_fML, c[i+1][j-1-i-1] + E_MLstem(tt, S1[i], S1[j], P) + 2*P->MLbase);
                  break;
      }

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

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

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

      /* coaxial stacking */
      if (dangle_model==3) {
        /* additional ML decomposition as two coaxially stacked helices */
        for (decomp = INF, k = i+1+turn; k <= j-2-turn; k++) {
          type = ptype[i][k-i]; type = rtype[type];
          type_2 = ptype[k+1][j-k-1]; type_2 = rtype[type_2];
          if (type && type_2)
            decomp = MIN2(decomp,
                          c[i][k-i]+c[k+1][j-k-1]+P->stack[type][type_2]);
        }

        decomp += 2*P->MLintern[1];          /* no TermAU penalty if coax stack */
#if 0
        /* This is needed for Y shaped ML loops with coax stacking of
           interior pairts, but backtracking will fail if activated */
        DMLi[j-i] = MIN2(DMLi[j-i], decomp);
        DMLi[j-i] = MIN2(DMLi[j-i], DMLi[j-1-i]+P->MLbase);
        DMLi[j-i] = MIN2(DMLi[j-i], DMLi1[j-(i+1)]+P->MLbase);
        new_fML = MIN2(new_fML, DMLi[j-i]);
#endif
        new_fML = MIN2(new_fML, decomp);
      }
      fML[i][j-i] = Fmi[j-i] = new_fML;     /* substring energy */
    } /* for (j...) */

    /* calculate energies of 5' and 3' fragments */
    {
      static int do_backtrack = 0, prev_i=0;
      char *ss=NULL;
      double prevz = 0.;

      /* first case: i stays unpaired */
      f3[i] = f3[i+1];

      /* next all cases where i is paired */
      switch(dangle_model){
        /* dont use dangling end and mismatch contributions at all */
        case 0:   for(j=i+turn+1; j<length && j<=i+maxdist; j++){
                    type = ptype[i][j-i];

                    if(with_gquad){
                      f3[i] = MIN2(f3[i], f3[j+1] + ggg[i][j-i]);
                    }

                    if(type)
                      f3[i] = MIN2(f3[i], f3[j+1] + c[i][j-i] + E_ExtLoop(type, -1, -1, P));
                  }
                  if(length<=i+maxdist){
                    j=length;

                    if(with_gquad){
                      f3[i] = MIN2(f3[i], ggg[i][j-i]);
                    }

                    type = ptype[i][j-i];
                    if(type)
                      f3[i] = MIN2(f3[i], c[i][j-i] + E_ExtLoop(type, -1, -1, P));
                  }
                  break;
        /* always use dangle_model on both sides */
        case 2:   for(j=i+turn+1; j<length && j<=i+maxdist; j++){
                    type = ptype[i][j-i];

                    if(with_gquad){
                      if(ggg[i][j-i] != INF)
                        f3[i] = MIN2(f3[i], f3[j+1] + ggg[i][j-i]);
                    }

                    if(type)
                      f3[i] = MIN2(f3[i], f3[j+1] + c[i][j-i] + E_ExtLoop(type, (i>1) ? S1[i-1] : -1, S1[j+1], P));
                  }
                  if(length<=i+maxdist){
                    j=length;

                    if(with_gquad){
                      f3[i] = MIN2(f3[i], ggg[i][j-i]);
                    }

                    type = ptype[i][j-i];
                    if(type)
                      f3[i] = MIN2(f3[i], c[i][j-i] + E_ExtLoop(type, (i>1) ? S1[i-1] : -1, -1, P));
                  }
                  break;
        /* normal dangle_model, aka dangle_model = 1 */
        default:  for(j=i+turn+1; j<length && j<=i+maxdist; j++){
                    type = ptype[i][j-i];

                    if(with_gquad){
                      f3[i] = MIN2(f3[i], f3[j+1] + ggg[i][j-i]);
                    }

                    if(type){
                      f3[i] = MIN2(f3[i], f3[j+1] + c[i][j-i] + E_ExtLoop(type, -1, -1, P));
                      f3[i] = MIN2(f3[i], ((j+2<=length) ? f3[j+2] : 0) + c[i][j-i] + E_ExtLoop(type, -1, S1[j+1], P));
                    }
                    type = ptype[i+1][j-i-1];
                    if(type){
                      f3[i] = MIN2(f3[i], f3[j+1] + c[i+1][j-i-1] + E_ExtLoop(type, S1[i], -1, P));
                      f3[i] = MIN2(f3[i], ((j + 1 < length) ? f3[j+2] : 0) + c[i+1][j-i-1] + E_ExtLoop(type, S1[i], S1[j+1], P));
                    }
                  }
                  if(length<=i+maxdist){
                    j     = length;

                    if(with_gquad){
                      f3[i] = MIN2(f3[i], ggg[i][j-i]);
                    }

                    type  = ptype[i][j-i];
                    if(type)
                      f3[i] = MIN2(f3[i], c[i][j-i] + E_ExtLoop(type, -1, -1, P));
                    type  = ptype[i+1][j-i-1];
                    if(type)
                      f3[i] = MIN2(f3[i], c[i+1][j-i-1] + E_ExtLoop(type, S1[i], -1, P));
                  }
                  break;
      } /* switch(dangle_model)... */

      /* backtrack partial structure */
      if (f3[i] < f3[i+1]){
        do_backtrack=1;
      }
      else if (do_backtrack) {
        int pairpartner; /*i+1?? is paired with pairpartner*/
        int cc;
        int traced2=0;
        fij = f3[i+1];
        lind=i+1;
        /*start "short" backtrack*/

        /*get paired base*/
        while(fij==f3[lind+1])
          lind++;

        /*get pairpartner*/
        for (pairpartner = lind + turn; pairpartner <= lind + maxdist; pairpartner++){
          type = ptype[lind][pairpartner-lind];
          switch(dangle_model){
            case 0:   if(type){
                        cc = c[lind][pairpartner-lind] + E_ExtLoop(type, -1, -1, P);
                        if(fij == cc + f3[pairpartner + 1])
                          traced2 = 1;
                      }
                      else if(with_gquad) {
                        cc = ggg[lind][pairpartner-lind];
                        if(fij == cc + f3[pairpartner + 1])
                          traced2 = 1;
                      }

                      break;
            case 2:   if(type){
                        cc = c[lind][pairpartner-lind] + E_ExtLoop(type, (lind > 1) ? S1[lind-1] : -1, (pairpartner < length) ? S1[pairpartner+1] : -1, P);
                        if(fij == cc + f3[pairpartner + 1])
                          traced2 = 1;
                      }
                      else if(with_gquad){
                        cc = ggg[lind][pairpartner-lind];
                        if(fij == cc + f3[pairpartner + 1])
                          traced2 = 1;
                      }

                      break;
            default:  if(type){
                        cc = c[lind][pairpartner-lind] + E_ExtLoop(type, -1, -1, P);
                        if(fij == cc + f3[pairpartner + 1]){
                          traced2 = 1;
                          break;
                        }
                        else if(pairpartner < length){
                          cc = c[lind][pairpartner-lind] + E_ExtLoop(type, -1, S1[pairpartner+1], P);
                          if(fij == cc + f3[pairpartner + 2]){
                            traced2 = 1;
                            break;
                          }
                        }
                      }
                      else if(with_gquad){
                        cc = ggg[lind][pairpartner-lind];
                        if(fij == cc + f3[pairpartner + 1])
                          traced2 = 1;
                      }

                      type = ptype[lind+1][pairpartner-lind-1];
                      if(type){
                        cc = c[lind+1][pairpartner-(lind+1)] + E_ExtLoop(type, S1[lind], -1, P);
                        if(fij == cc + f3[pairpartner+1]){
                          traced2 = 1;
                          break;
                        }
                        else if(pairpartner < length){
                          cc = c[lind+1][pairpartner-(lind+1)] + E_ExtLoop(type, S1[lind], S1[pairpartner+1], P);
                          if(fij == cc + f3[pairpartner+2])
                            traced2 = 1;
                        }
                      }
                      break;
          }
          if(traced2) break;
        }
        if (!traced2) vrna_message_error("backtrack failed in short backtrack 1");
        if (zsc){
#ifdef USE_SVM
          int info_avg;
          double average_free_energy;
          double sd_free_energy;
          double my_z;
          int *AUGC = get_seq_composition(S, lind-1, MIN2((pairpartner+1),length), length);
          /*\svm*/
          average_free_energy = avg_regression(AUGC[0], AUGC[1], AUGC[2], AUGC[3], AUGC[4], avg_model, &info_avg);
          if (info_avg == 0)  {
            double difference;
            double min_sd = minimal_sd(AUGC[0],AUGC[1],AUGC[2],AUGC[3],AUGC[4]);
            difference=(fij-f3[pairpartner+1])/100.-average_free_energy;
            if ( difference - ( min_z * min_sd ) <= 0.0001 ) {
              sd_free_energy = sd_regression(AUGC[0],AUGC[1],AUGC[2],AUGC[3],AUGC[4],sd_model);
              my_z=difference/sd_free_energy;
              if (my_z<=min_z){
                ss =  backtrack(vc, lind, pairpartner+1);
                if (prev) {
                  if ((i+strlen(ss)<prev_i+strlen(prev)) ||
                      strncmp(ss+prev_i-i,prev,strlen(prev))) { /* ss does not contain prev */
                    if (dangle_model==2)
                      fprintf(output, ".%s (%6.2f) %4d z= %.3f\n", prev, (f3[prev_i]-f3[prev_i+strlen(prev)-1])/100., prev_i-1, prevz);
                    else
                      fprintf(output, "%s (%6.2f) %4d z=%.3f\n ", prev, (f3[prev_i]-f3[prev_i+strlen(prev)])/100., prev_i, prevz);
                  }
                  free(prev);
                }
                prev=ss; prev_i = lind; prevz=my_z;
              }
            }

          }
          free(AUGC);
          do_backtrack=0;
#endif
        }
        else {
          /* original code for Lfold*/
          ss =  backtrack(vc, lind , pairpartner+1);
          if (prev) {
            if ((i+strlen(ss)<prev_i+strlen(prev)) || strncmp(ss+prev_i-i,prev,strlen(prev))){
              /* ss does not contain prev */
              if (dangle_model==2){
                fprintf(output, ".%s (%6.2f) %4d\n", prev, (f3[prev_i]-f3[prev_i+strlen(prev)-1])/100., prev_i-1);
              } else
                fprintf(output, "%s (%6.2f) %4d\n", prev, (f3[prev_i]-f3[prev_i+strlen(prev)])/100., prev_i);
            }
            free(prev);
          }
          prev=ss;
          prev_i = lind;
          do_backtrack=0;
        }
      }
      if (i==1) {
        if (prev) {
          if(zsc) {
            if (dangle_model==2)
              fprintf(output, ".%s (%6.2f) %4d z= %.2f\n", prev, (f3[prev_i]-f3[prev_i+strlen(prev)-1])/100., prev_i-1, prevz);
           else
              fprintf(output, "%s (%6.2f) %4dz= %.2f \n", prev, (f3[prev_i]-f3[prev_i+strlen(prev)])/100., prev_i, prevz);
          }
          else {
            if (dangle_model==2)
              fprintf(output, ".%s (%6.2f) %4d\n", prev, (f3[prev_i]-f3[prev_i+strlen(prev)-1])/100., prev_i-1);
            else
              fprintf(output, "%s (%6.2f) %4d\n", prev, (f3[prev_i]-f3[prev_i+strlen(prev)])/100., prev_i);
          }
          free(prev); prev=NULL;
        } else if ((f3[i]<0) && (!zsc)) do_backtrack=1;

        if (do_backtrack) {
          int pairpartner; /*i+1?? is paired with pairpartner*/
          int cc;
          double average_free_energy;
          double sd_free_energy;
          int info_avg;
          double my_z;
          int traced2 = 0;
          fij = f3[i];
          lind=i;
          while(fij==f3[lind+1]) lind++;
          /*get pairpartner*/
          for(pairpartner = lind + turn; pairpartner <= lind + maxdist; pairpartner++){
            type = ptype[lind][pairpartner-lind];
            switch(dangle_model){
              case 0:   if(type){
                          cc = c[lind][pairpartner-lind] + E_ExtLoop(type, -1, -1, P);
                          if(fij == cc + f3[pairpartner + 1])
                            traced2 = 1;
                        }
                        else if(with_gquad){
                          cc = ggg[lind][pairpartner-lind];
                          if(fij == cc + f3[pairpartner + 1])
                            traced2 = 1;
                        }

                        break;
              case 2:   if(type){
                          cc = c[lind][pairpartner-lind] + E_ExtLoop(type, (lind > 1) ? S1[lind-1] : -1, (pairpartner < length) ? S1[pairpartner+1] : -1, P);
                          if(fij == cc + f3[pairpartner + 1])
                            traced2 = 1;
                        }
                        else if(with_gquad){
                          cc = ggg[lind][pairpartner-lind];
                          if(fij == cc + f3[pairpartner + 1])
                            traced2 = 1;
                        }

                        break;
              default:  if(type){
                          cc = c[lind][pairpartner-lind] + E_ExtLoop(type, -1, -1, P);
                          if(fij == cc + f3[pairpartner + 1]){
                            traced2 = 1;
                            break;
                          }
                          else if(pairpartner < length){
                            cc = c[lind][pairpartner-lind] + E_ExtLoop(type, -1, S1[pairpartner + 1], P);
                            if(fij == cc + f3[pairpartner + 1]){
                              traced2 = 1;
                              break;
                            }
                          }
                        }
                        else if(with_gquad){
                          cc = ggg[lind][pairpartner-lind];
                          if(fij == cc + f3[pairpartner + 1])
                            traced2 = 1;
                        }

                        type = ptype[lind+1][pairpartner-lind-1];
                        if(type){
                          cc = c[lind+1][pairpartner-(lind+1)] + E_ExtLoop(type, S1[lind], -1, P);
                          if(fij == cc + f3[pairpartner+1]){
                            traced2 = 1;
                            break;
                          }
                          else if (pairpartner < length){
                            cc = c[lind+1][pairpartner-(lind+1)] + E_ExtLoop(type, S1[lind], S1[pairpartner+1], P);
                            if(fij == cc + f3[pairpartner + 2]){
                              traced2 =1;
                              break;
                            }
                          }
                        }
            }
            if(traced2) break;
          }
          if (!traced2) vrna_message_error("backtrack failed in short backtrack 2");

          if(zsc){
#ifdef USE_SVM
            int *AUGC = get_seq_composition(S, lind-1, MIN2((pairpartner+1),length), length);
            average_free_energy = avg_regression(AUGC[0],AUGC[1],AUGC[2],AUGC[3],AUGC[4],avg_model,&info_avg);
            if (info_avg == 0)  {
              double difference;
              double min_sd = minimal_sd(AUGC[0],AUGC[1],AUGC[2],AUGC[3],AUGC[4]);
              difference=(fij-f3[pairpartner+1])/100.-average_free_energy;
              if ( difference - ( min_z * min_sd ) <= 0.0001 ) {
                sd_free_energy = sd_regression(AUGC[0],AUGC[1],AUGC[2],AUGC[3],AUGC[4],sd_model);
                my_z=difference/sd_free_energy;
                if (my_z<=min_z){
                  ss =  backtrack(vc, lind , pairpartner+1);
                  fprintf(output, "%s (%6.2f) %4d z= %.2f\n", ss, (f3[lind]-f3[lind+strlen(ss)-1])/100., lind, my_z);
                }
              }
            }
            free(AUGC);
#endif
          }
          else {
            ss =  backtrack(vc, lind , pairpartner+1);
            if (dangle_model==2)
              fprintf(output, "%s (%6.2f) %4d\n", ss, (f3[lind]-f3[lind+strlen(ss)-1])/100., 1);
            else
              fprintf(output, "%s (%6.2f) %4d\n", ss, (f3[lind]-f3[lind+strlen(ss)])/100., 1);
            free(ss);
          }
        }
        do_backtrack=0;
      }
    }
    {
      int ii, *FF; /* rotate the auxilliary arrays */

      /* check for values close to integer underflow */
      if(INT_CLOSE_TO_UNDERFLOW(f3[i])){
        /* correct f3 free energies and increase underflow counter */
        int cnt, cnt2;
        for(cnt=i; cnt <= length && cnt <= lind + maxdist + 2; cnt++) {
          f3[cnt] -= UNDERFLOW_CORRECTION;
        }
        (*underflow)++;
      }

      FF = DMLi2; DMLi2 = DMLi1; DMLi1 = DMLi; DMLi = FF;
      FF = cc1; cc1=cc; cc=FF;
      for(j = 0; j < maxdist + 5; j++){
        cc[j] = Fmi[j] = DMLi[j] = INF;
      }

      /*
        rotate the DP matrices
        NOTE: here we rotate them only locally, i.e. their
        actual configuration within vc remains intact
      */
      if( i + maxdist + 4 <= length ){
        c[i - 1]                = c[i + maxdist + 4];
        c[i + maxdist + 4]      = NULL;
        fML[i - 1]              = fML[i + maxdist + 4];
        fML[i + maxdist + 4]    = NULL;
        ptype[i - 1]            = ptype[i + maxdist + 4];
        ptype[i + maxdist + 4]  = NULL;
        if( i > 1 ){
          make_ptypes(vc, i - 1);
          if(with_gquad){
            vrna_gquad_mx_local_update(vc, i - 1);
            ggg = vc->matrices->ggg_local;
          }
        }
        for(ii = 0; ii < maxdist + 5; ii++){
          c[i - 1][ii]    = INF;
          fML[i - 1][ii]  = INF;
        }
      }

    }
  }

  free(cc);
  free(cc1);
  free(Fmi);
  free(DMLi);
  free(DMLi1);
  free(DMLi2);

  return f3[1];
}

PRIVATE char *
backtrack(vrna_fold_compound_t *vc,
          int start,
          int maxdist){

  /*------------------------------------------------------------------
    trace back through the "c", "f3" and "fML" arrays to get the
    base pairing list. No search for equivalent structures is done.
    This is fast, since only few structure elements are recalculated.
    ------------------------------------------------------------------*/
  sect          sector[MAXSECTORS];   /* backtracking sectors */
  int           i, j, k, length, energy, new, no_close, type, type_2, tt, s=0;
  int           with_gquad, bt_type, turn, dangle_model, noLP, noGUclosure, *rtype;
  int           **c, **fML, *f3, **ggg;
  char          *string, *structure, **ptype;
  short         *S, *S1;
  vrna_param_t  *P;
  vrna_md_t     *md;

  string        = vc->sequence;
  length        = vc->length;
  S             = vc->sequence_encoding2;
  S1            = vc->sequence_encoding;
  ptype         = vc->ptype_local;
  P             = vc->params;
  md            = &(P->model_details);
  dangle_model  = md->dangles;
  noLP          = md->noLP;
  noGUclosure   = md->noGUclosure;
  with_gquad    = md->gquad;
  bt_type       = md->backtrack_type;
  turn          = md->min_loop_size;
  rtype         = &(md->rtype[0]);

  c       = vc->matrices->c_local;
  fML     = vc->matrices->fML_local;
  f3      = vc->matrices->f3_local;
  ggg     = vc->matrices->ggg_local;

  /* length = strlen(string); */
  sector[++s].i = start;
  sector[s].j   = MIN2(length, maxdist+1);
  sector[s].ml  = (bt_type=='M') ? 1 : ((bt_type=='C')?2:0);

  structure = (char *) vrna_alloc((MIN2(length-start, maxdist)+3)*sizeof(char));
  for (i=0; i<=MIN2(length-start, maxdist); i++) structure[i] = '-';

  while (s>0) {
    int ml, fij, cij, traced, i1, j1, mm, mm5, mm3, mm53, p, q, 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) {
      structure[i-start] = '(';
      structure[j-start] = ')';
      goto repeat1;
    }

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

    fij = (ml)? fML[i][j-i] : f3[i];

    if (ml == 0) { /* backtrack in f3 */

      if (fij == f3[i+1]) {
        sector[++s].i = i+1;
        sector[s].j   = j;
        sector[s].ml  = ml;
        continue;
      }
      /* i or i+1 is paired. Find pairing partner */
      switch(dangle_model){
        case 0:   for(traced = 0, k=j; k>i+turn; k--){

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

                    jj    = k+1;
                    type  = ptype[i][k-i];
                    if(type)
                      if(fij == c[i][k-i] + E_ExtLoop(type, -1, -1, P) + f3[k+1]){
                        traced = i;
                        break;
                      }
                  }
                  break;
        case 2:   for(traced = 0, k=j; k>i+turn; k--){

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

                    jj    = k+1;
                    type  = ptype[i][k-i];
                    if(type)
                      if(fij == c[i][k-i] + E_ExtLoop(type, (i>1) ? S1[i-1] : -1, (k<length) ? S1[k+1] : -1, P) + f3[k+1]){
                        traced = i;
                        break;
                      }
                  }
                  break;
        default:  for(traced = 0,k=j; k>i+turn; k--){

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

                    jj = k+1;
                    type = ptype[i+1][k-(i+1)];
                    if(type){
                      if(fij == c[i+1][k-(i+1)] + E_ExtLoop(type, S1[i], -1, P) + f3[k+1]){
                        traced=i+1;
                      }
                      if(k < length){
                        if(fij == c[i+1][k-(i+1)] + E_ExtLoop(type, S1[i], S1[k+1], P) + f3[k+2]){
                          traced  = i+1;
                          jj      = k+2;
                        }
                      }
                    }
                    type = ptype[i][k-i];
                    if(type){
                      if(fij == c[i][k-i] + E_ExtLoop(type, -1, -1, P) + f3[k+1]){
                        traced = i;
                      }
                      if(k<length){
                        if(fij == c[i][k-i] + E_ExtLoop(type, -1, S1[k+1], P) + f3[k+2]){
                          traced  = i;
                          jj      = k+2;
                        }
                      }
                    }
                    if(traced) break;
                  }
                  break;
      } /* switch(dangle_model)...*/

      if (!traced) vrna_message_error("backtrack failed in f3");
      if (j==length) { /* backtrack only one component, unless j==length */
        sector[++s].i = jj;
        sector[s].j   = j;
        sector[s].ml  = ml;
      }
      i=traced; j=k;

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

      structure[i-start] = '('; structure[j-start] = ')';
      if (((jj==j+2) || (dangle_model==2)) && (j < length)) structure[j+1-start] = '.';
      goto repeat1;
    }
    else { /* trace back in fML array */
      if (fML[i][j-1-i]+P->MLbase == fij) {  /* 3' end is unpaired */
        sector[++s].i = i;
        sector[s].j   = j-1;
        sector[s].ml  = ml;
        continue;
      }
      if (fML[i+1][j-(i+1)]+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[i][j-i] + E_MLstem(0, -1, -1, P)){
          /* go to backtracing of quadruplex */
          goto repeat_gquad;
        }
      }

      switch(dangle_model){
        case 0:   tt = ptype[i][j-i];
                  if(fij == c[i][j-i] + E_MLstem(tt, -1, -1, P)){
                    structure[i-start] = '(';
                    structure[j-start] = ')';
                    goto repeat1;
                  }
                  break;
        case 2:   tt = ptype[i][j-i];
                  if(fij == c[i][j-i] + E_MLstem(tt, (i>1) ? S1[i-1] : -1, (j < length) ? S1[j+1] : -1, P)){
                    structure[i-start] = '(';
                    structure[j-start] = ')';
                    goto repeat1;
                  }
                  break;
        default:  tt = ptype[i][j-i];
                  if(fij == c[i][j-i] + E_MLstem(tt, -1, -1, P)){
                    structure[i-start] = '(';
                    structure[j-start] = ')';
                    goto repeat1;
                  }
                  tt = ptype[i+1][j-(i+1)];
                  if(fij == c[i+1][j-(i+1)] + E_MLstem(tt, S1[i], -1, P) + P->MLbase){
                    structure[++i-start] = '(';
                    structure[j-start] = ')';
                    goto repeat1;
                  }
                  tt = ptype[i][j-1-i];
                  if(fij == c[i][j-1-i] + E_MLstem(tt, -1, S1[j], P) + P->MLbase){
                    structure[i-start] = '(';
                    structure[--j-start] = ')';
                    goto repeat1;
                  }
                  tt = ptype[i+1][j-1-(i+1)];
                  if(fij == c[i+1][j-1-(i+1)] + E_MLstem(tt, S1[i], S1[j], P) + 2*P->MLbase){
                    structure[++i-start] = '(';
                    structure[--j-start] = ')';
                    goto repeat1;
                  }
                  break;
      } /* switch(dangle_model)... */

      /* modular decomposition */
      for (k = i+1+turn; k <= j-2-turn; k++)
        if (fij == (fML[i][k-i]+fML[k+1][j-(k+1)]))
          break;

      if ((dangle_model==3)&&(k>j-2-turn)) { /* must be coax stack */
        ml = 2;
        for (k = i+1+turn; k <= j-2-turn; k++) {
          type = ptype[i][k-i];  type= rtype[type];
          type_2 = ptype[k+1][j-(k+1)]; type_2= rtype[type_2];
          if (type && type_2)
            if (fij == c[i][k-i]+c[k+1][j-(k+1)]+P->stack[type][type_2]+
                       2*P->MLintern[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) vrna_message_error("backtrack failed in fML");
      continue;
    }

  repeat1:

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

    type = ptype[i][j-i];


    if (noLP)
      if (cij == c[i][j-i]) {
        /* (i.j) closes canonical structures, thus
           (i+1.j-1) must be a pair                */
        type_2 = ptype[i+1][j-1-(i+1)]; type_2 = rtype[type_2];
        cij -= P->stack[type][type_2];
        structure[i+1-start] = '('; structure[j-1-start] = ')';
        i++; j--;
        canonical=0;
        goto repeat1;
      }
    canonical = 1;


    no_close = (((type==3)||(type==4))&&noGUclosure);
    if (no_close) {
      if (cij == FORBIDDEN) continue;
    } else
      if (cij == E_Hairpin(j-i-1, type, S1[i+1], S1[j-1],string+i-1, P))
        continue;

    for (p = i+1; p <= MIN2(j-2-turn,i+MAXLOOP+1); p++) {
      int minq;
      minq = j-i+p-MAXLOOP-2;
      if (minq<p+1+turn) minq = p+1+turn;
      for (q = j-1; q >= minq; q--) {

        type_2 = ptype[p][q-p];
        if (type_2==0) continue;
        type_2 = rtype[type_2];
        if (noGUclosure)
          if (no_close||(type_2==3)||(type_2==4))
            if ((p>i+1)||(q<j-1)) continue;  /* continue unless stack */

        /* energy = oldLoopEnergy(i, j, p, q, type, type_2); */
        energy = E_IntLoop(p-i-1, j-q-1, type, type_2, S1[i+1], S1[j-1], S1[p-1], S1[q+1],P);

        new = energy+c[p][q-p];
        traced = (cij == new);
        if (traced) {
          structure[p-start] = '(';
          structure[q-start] = ')';
          i = p, j = q;
          goto repeat1;
        }
      }
    }

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

    /* (i.j) must close a multi-loop */
    tt = rtype[type];
    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...
      */
      if(backtrack_GQuad_IntLoop_L(cij, i, j, type, S, ggg, maxdist, &p, &q, P)){
        i = p; j = q;
        goto repeat_gquad;
      }
    }

    sector[s+1].ml  = sector[s+2].ml = 1;

    switch(dangle_model){
      case 0:   mm = P->MLclosing + E_MLstem(tt, -1, -1, P);
                for(k = i+2+turn; k < j-2-turn; k++){
                  if(cij == fML[i+1][k-(i+1)] + fML[k+1][j-1-(k+1)] + mm)
                    break;
                }
                break;
      case 2:   mm = P->MLclosing + E_MLstem(tt, S1[j-1], S1[i+1], P);
                for(k = i+2+turn; k < j-2-turn; k++){
                  if(cij == fML[i+1][k-(i+1)] + fML[k+1][j-1-(k+1)] + mm)
                    break;
                }
                break;
      default:  mm    = P->MLclosing + E_MLstem(tt, -1, -1, P);
                mm5   = P->MLclosing + E_MLstem(tt, S1[j-1], -1, P) + P->MLbase;
                mm3   = P->MLclosing + E_MLstem(tt, -1, S1[i+1], P) + P->MLbase;
                mm53  = P->MLclosing + E_MLstem(tt, S1[j-1], S1[i+1], P) + 2*P->MLbase;
                for(k = i+2+turn; k < j-2-turn; k++){
                  if(cij == fML[i+1][k-(i+1)] + fML[k+1][j-1-(k+1)] + mm)
                    break;
                  else if(cij == fML[i+2][k-(i+2)] + fML[k+1][j-1-(k+1)] + mm3){
                    i1 = i+2;
                    break;
                  }
                  else if(cij == fML[i+1][k-(i+1)] + fML[k+1][j-2-(k+1)] + mm5){
                    j1 = j-2;
                    break;
                  }
                  else if(cij == fML[i+2][k-(i+2)] + fML[k+1][j-2-(k+1)] + mm53){
                    i1 = i+2;
                    j1 = j-2;
                    break;
                  }
                  /* coaxial stacking of (i.j) with (i+1.k) or (k.j-1) */
                  /* use MLintern[1] since coax stacked pairs don't get TerminalAU */
                  if (dangle_model==3) {
                    int en;
                    type_2 = ptype[i+1][k-(i+1)]; type_2 = rtype[type_2];
                    if (type_2) {
                      en = c[i+1][k-(i+1)]+P->stack[type][type_2]+fML[k+1][j-1-(k+1)];
                      if (cij == en+2*P->MLintern[1]+P->MLclosing) {
                        ml = 2;
                        sector[s+1].ml  = 2;
                        break;
                      }
                    }
                    type_2 = ptype[k+1][j-1-(k+1)]; type_2 = rtype[type_2];
                    if (type_2) {
                      en = c[k+1][j-1-(k+1)]+P->stack[type][type_2]+fML[i+1][k-(i+1)];
                      if (cij == en+2*P->MLintern[1]+P->MLclosing) {
                        sector[s+2].ml = 2;
                        break;
                      }
                    }
                  }
                }
                break;
    } /* switch(dangle_model)... */

    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 {
#if 0
      /* Y shaped ML loops fon't work yet */
      if (dangle_model==3) {
        /* (i,j) must close a Y shaped ML loop with coax stacking */
        if (cij ==  fML[i+1][j-2-(i+2)] + mm + d3 + d5 + P->MLbase + P->MLbase) {
          i1 = i+2;
          j1 = j-2;
        } else if (cij ==  fML[i+1][j-2-(i+1)] + mm + d5 + P->MLbase)
          j1 = j-2;
        else if (cij ==  fML[i+2][j-1-(i+2)] + mm + d3 + P->MLbase)
          i1 = i+2;
        else /* last chance */
          if (cij != fML[i+1][j-1-(i+1)] + mm + P->MLbase)
            fprintf(stderr,  "backtracking failed in repeat");
        /* if we arrive here we can express cij via fML[i1,j1]+dangle_model */
        sector[++s].i = i1;
        sector[s].j   = j1;
      }
      else
#endif
        vrna_message_error("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 l[3], L, a;
      L = -1;

      get_gquad_pattern_mfe(S, i, j, P, &L, l);
      if(L != -1){
        /* fill the G's of the quadruplex into the structure string */
        for(a=0;a<L;a++){
          structure[i+a-start] = '+';
          structure[i+L+l[0]+a-start] = '+';
          structure[i+L+l[0]+L+l[1]+a-start] = '+';
          structure[i+L+l[0]+L+l[1]+L+l[2]+a-start] = '+';
        }
        goto repeat_gquad_exit;
      }
      vrna_message_error("backtracking failed in repeat_gquad");
    }
  repeat_gquad_exit:
    __asm("nop");

  }

  for (i=strlen(structure)-1; i>0 && structure[i] == '-'; i--)
    structure[i] = '\0';
  for (;i>=0; i--)
   if (structure[i]=='-') structure[i]='.';

  return structure;
}

PRIVATE void
make_ptypes(vrna_fold_compound_t *vc, int i){

  int       j, k, type, n, maxdist, turn, noLP;
  short     *S;
  char      **ptype;
  vrna_md_t *md;

  n       = (int)vc->length;
  S       = vc->sequence_encoding2;
  ptype   = vc->ptype_local;
  maxdist = vc->window_size;
  md      = &(vc->params->model_details);
  turn    = md->min_loop_size;
  noLP    = md->noLP;

  for(k = turn + 1; k < maxdist; k++){
    j = i + k;
    if (j > n)
      break;
    type = md->pair[S[i]][S[j]];

    if(noLP && type){
      if(!ptype[i + 1][j - 1 - i - 1])
        if(j == n || i == 1 || (!md->pair[S[i - 1]][S[j + 1]]))
          type = 0;
    }
    ptype[i][j - i] = type;
  }
}

/*###########################################*/
/*# deprecated functions below              #*/
/*###########################################*/

#ifdef  VRNA_BACKWARD_COMPAT

PUBLIC float Lfold( const char *string,
                    char *structure,
                    int window_size){

  float               energy;
  vrna_fold_compound_t  *vc;
  vrna_md_t           md;

  set_model_details(&md);

  md.window_size = window_size;
  md.max_bp_span = window_size;

  vc  = vrna_fold_compound(string, &md, VRNA_OPTION_WINDOW);

  energy = wrap_Lfold(vc, 0, 0.0, NULL);

  vrna_fold_compound_free(vc);

  return energy;
}

PUBLIC float
Lfoldz( const char *string,
        char *structure,
        int window_size,
        int zsc,
        double min_z){

  float               energy;
  vrna_fold_compound_t  *vc;
  vrna_md_t           md;

  set_model_details(&md);

  md.window_size = window_size;
  md.max_bp_span = window_size;

  vc  = vrna_fold_compound(string, &md, VRNA_OPTION_WINDOW);

#ifndef USE_SVM
  zsc = 0;  /* deactivate z-scoring if no compiled-in svm support is available */
#endif

  energy = wrap_Lfold(vc, zsc, min_z, NULL);

  vrna_fold_compound_free(vc);

  return energy;
}

#endif