/*
  $Log: subopt.c,v $
  Revision 2.0  2010/12/06 20:04:20  ronny
  repaired subopt for cofolding

  Revision 1.24  2008/11/01 21:10:20  ivo
  avoid rounding errors when computing DoS

  Revision 1.23  2008/03/31 15:06:49  ivo
  Add cofolding support in subopt

  Revision 1.22  2008/02/23 09:42:35  ivo
  fix circular folding bugs with dangles that cross the origin

  Revision 1.21  2008/01/08 15:08:51  ivo
  circular fold would fail for open chain

  Revision 1.20  2008/01/08 14:08:20  ivo
  add an option to compute the density of state

  Revision 1.19  2007/12/05 13:04:04  ivo
  add various circfold variants from Ronny

  Revision 1.18  2003/10/06 08:56:45  ivo
  use P->TerminalAU

  Revision 1.17  2003/08/26 09:26:08  ivo
  don't modify print_energy in subopt(); use doubles instead of floats

  Revision 1.16  2001/10/01 13:50:00  ivo
  sorted -> subopt_sorted

  Revision 1.15  2001/09/17 10:30:42  ivo
  move scale_parameters() into params.c
  returns pointer to paramT structure

  Revision 1.14  2001/08/31 15:02:19  ivo
  Let subopt either write to file pointer or return a list of structures,
  so we can nicely integrate it into the library

  Revision 1.13  2001/04/05 07:35:08  ivo
  remove uneeded declaration of TETRA_ENERGY

  Revision 1.12  2000/10/10 08:53:20  ivo
  adapted for new Turner energy parameters
  supports all constraints that forbid pairs

  Revision 1.11  2000/04/08 15:56:18  ivo
  with noLonelyPairs=1 will produce no structures with isolated base pairs
  (Giegerich's canonical structures)

  Revision 1.10  1999/05/06 10:13:35  ivo
  recalculte energies before printing if logML is set
  + cosmetic changes

  Revision 1.9  1998/05/19 16:31:52  ivo
  added support for constrained folding

  Revision 1.8  1998/03/30 14:44:54  ivo
  cleanup of make_printout etc.

  Revision 1.7  1998/03/30 14:39:31  ivo
  replaced BasePairs list with structure string in STATE
  save memory by not storing (and sorting) structures
  modified for use with ViennaRNA-1.2.1

  Revision 1.6  1997/10/21 11:34:09  walter
  steve update

  Revision 1.1  1997/08/04 21:05:32  walter
  Initial revision

*/
/*
   suboptimal folding - Stefan Wuchty, Walter Fontana & Ivo Hofacker

                       Vienna RNA package
*/
#include <config.h>
#include <stdio.h>
#include <stdlib.h>
#include <ctype.h>
#include <string.h>
#include <math.h>
#include "fold.h"
#include "utils.h"
#include "energy_par.h"
#include "fold_vars.h"
#include "pair_mat.h"
#include "list.h"
#include "params.h"
#include "loop_energies.h"
#include "cofold.h"
#include "subopt.h"

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

#define true              1
#define false             0
#define SAME_STRAND(I,J)  (((I)>=cut_point)||((J)<cut_point))
#define NEW_NINIO         1         /* use new asymetry penalty */
#define STACK_BULGE1      1         /* stacking energies for bulges of size 1 */
#define MAXALPHA          20        /* maximal length of alphabet */

/*@unused@*/
PRIVATE char UNUSED rcsid[] = "$Id: subopt.c,v 1.24 2008/11/01 21:10:20 ivo Exp $";

/*
#################################
# GLOBAL VARIABLES              #
#################################
*/
PUBLIC  int     subopt_sorted=0;                           /* output sorted by energy */
PUBLIC  int     density_of_states[MAXDOS+1];
PUBLIC  double  print_energy = 9999; /* printing threshold for use with logML */

typedef struct {
    char *structure;
    LIST *Intervals;
    int partial_energy;
    int is_duplex;
    /* int best_energy;   */ /* best attainable energy */
} STATE;

/*
#################################
# PRIVATE VARIABLES             #
#################################
*/
PRIVATE int     turn;
PRIVATE LIST    *Stack = NULL;
PRIVATE int     nopush;
PRIVATE int     best_energy;          /* best_energy = remaining energy */
PRIVATE int     *f5 = NULL;           /* energy of 5 end */
PRIVATE int     *c = NULL;            /* energy array, given that i-j pair */
PRIVATE int     *fML = NULL;          /* multi-loop auxiliary energy array */
PRIVATE int     *fM1 = NULL;          /* another multi-loop auxiliary energy array */
PRIVATE int     *fc = NULL;           /*energy array, from i (j)  to cut*/
PRIVATE int     *indx = NULL;         /* index for moving in the triangle matrices c[] and f[] */
PRIVATE short   *S=NULL, *S1=NULL;
PRIVATE char    *ptype=NULL;
PRIVATE paramT  *P = NULL;
PRIVATE int     length;
PRIVATE int     minimal_energy;       /* minimum free energy */
PRIVATE int     element_energy;       /* internal energy of a structural element */
PRIVATE int     threshold;            /* minimal_energy + delta */
PRIVATE char    *sequence = NULL;
PRIVATE int     circular            = 0;
PRIVATE int     struct_constrained  = 0;
PRIVATE int     *fM2 = NULL;                 /* energies of M2 */
PRIVATE int     Fc, FcH, FcI, FcM;    /* parts of the exterior loop energies */

#ifdef _OPENMP

/* NOTE: all variables are assumed to be uninitialized if they are declared as threadprivate
         thus we have to initialize them before usage by a seperate function!
         OR: use copyin() in the PARALLEL directive!
*/
#pragma omp threadprivate(turn, Stack, nopush, best_energy, f5, c, fML, fM1, fc, indx, S, S1,\
                          ptype, P, length, minimal_energy, element_energy, threshold, sequence,\
                          fM2, Fc, FcH, FcI, FcM, circular, struct_constrained)

#endif

/*
#################################
# PRIVATE FUNCTION DECLARATIONS #
#################################
*/
PRIVATE void      make_pair(int i, int j, STATE *state);
PRIVATE INTERVAL  *make_interval (int i, int j, int ml);
/*@out@*/ PRIVATE STATE *make_state(/*@only@*/LIST *Intervals,
                                    /*@only@*/ /*@null@*/ char *structure,
                                    int partial_energy, int is_duplex);
PRIVATE STATE     *copy_state(STATE * state);
PRIVATE void      print_state(STATE * state);
PRIVATE void      UNUSED print_stack(LIST * list);
/*@only@*/ PRIVATE LIST *make_list(void);
PRIVATE void      push(LIST * list, /*@only@*/ void *data);
PRIVATE void      *pop(LIST * list);
PRIVATE int       best_attainable_energy(STATE * state);
PRIVATE void      scan_interval(int i, int j, int array_flag, STATE * state);
PRIVATE void      free_interval_node(/*@only@*/ INTERVAL * node);
PRIVATE void      free_state_node(/*@only@*/ STATE * node);
PRIVATE void      push_back(STATE * state);
PRIVATE char*     get_structure(STATE * state);
PRIVATE int       compare(const void *solution1, const void *solution2);
PRIVATE void      make_output(SOLUTION *SL, FILE *fp);
PRIVATE char      *costring(char *string);
PRIVATE void      repeat(int i, int j, STATE * state,
                  int part_energy, int temp_energy);

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



/*---------------------------------------------------------------------------*/
/*List routines--------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

PRIVATE void
make_pair(int i, int j, STATE *state)
{
  state->structure[i-1] = '(';
  state->structure[j-1] = ')';
}

/*---------------------------------------------------------------------------*/

PRIVATE INTERVAL *
make_interval(int i, int j, int array_flag)
{
  INTERVAL *interval;

  interval = lst_newnode(sizeof(INTERVAL));
  interval->i = i;
  interval->j = j;
  interval->array_flag = array_flag;
  return interval;
}

/*---------------------------------------------------------------------------*/

PRIVATE void
free_interval_node(INTERVAL * node)
{
  lst_freenode(node);
}

/*---------------------------------------------------------------------------*/

PRIVATE void
free_state_node(STATE * node)
{
  free(node->structure);
  if (node->Intervals)
    lst_kill(node->Intervals, lst_freenode);
  lst_freenode(node);
}

/*---------------------------------------------------------------------------*/

PRIVATE STATE *
make_state(LIST * Intervals,
           char *structure,
           int partial_energy,
           int is_duplex)
{
  STATE *state;

  state = lst_newnode(sizeof(STATE));

  if (Intervals)
    state->Intervals = Intervals;
  else
    state->Intervals = lst_init();
  if (structure)
    state->structure = structure;
  else {
    int i;
    state->structure = (char *) space(length+1);
    for (i=0; i<length; i++)
      state->structure[i] = '.';
  }

  state->partial_energy = partial_energy;

  return state;
}

/*---------------------------------------------------------------------------*/

PRIVATE STATE *
copy_state(STATE * state)
{
  STATE *new_state;
  void *after;
  INTERVAL *new_interval, *next;

  new_state = lst_newnode(sizeof(STATE));
  new_state->Intervals = lst_init();
  new_state->partial_energy = state->partial_energy;
  /* new_state->best_energy = state->best_energy; */

  if (state->Intervals->count) {
    after = LST_HEAD(new_state->Intervals);
    for ( next = lst_first(state->Intervals); next; next = lst_next(next))
      {
        new_interval = lst_newnode(sizeof(INTERVAL));
        *new_interval = *next;
        lst_insertafter(new_state->Intervals, new_interval, after);
        after = new_interval;
      }
  }
  new_state->structure = strdup(state->structure);
  if (!new_state->structure) nrerror("out of memory");
  return new_state;
}

/*---------------------------------------------------------------------------*/

/*@unused @*/ PRIVATE void
print_state(STATE * state)
{
  INTERVAL *next;

  if (state->Intervals->count)
    {
      printf("%d intervals:\n", state->Intervals->count);
      for (next = lst_first(state->Intervals); next; next = lst_next(next))
        {
          printf("[%d,%d],%d ", next->i, next->j, next->array_flag);
        }
      printf("\n");
    }
  printf("partial structure: %s\n", state->structure);
  printf("\n");
  printf(" partial_energy: %d\n", state->partial_energy);
  /* printf(" best_energy: %d\n", state->best_energy); */
  (void) fflush(stdout);
}

/*---------------------------------------------------------------------------*/

/*@unused @*/ PRIVATE void
print_stack(LIST * list)
{
  void *rec;

  printf("================\n");
  printf("%d states\n", list->count);
  for (rec = lst_first(list); rec; rec = lst_next(rec))
    {
      printf("state-----------\n");
      print_state(rec);
    }
  printf("================\n");
}

/*---------------------------------------------------------------------------*/

PRIVATE LIST *
make_list(void)
{
  return lst_init();
}

/*---------------------------------------------------------------------------*/

PRIVATE void
push(LIST * list, void *data)
{
  nopush = false;
  lst_insertafter(list, data, LST_HEAD(list));
}

/* PRIVATE void */
/* push_stack(STATE *state) { */ /* keep the stack sorted by energy */
/*   STATE *after, *next; */
/*   nopush = false; */
/*   next = after = LST_HEAD(Stack); */
/*   while ( next = lst_next(next)) { */
/*     if ( next->best_energy >= state->best_energy ) break; */
/*     after = next; */
/*   } */
/*   lst_insertafter(Stack, state, after); */
/* } */

/*---------------------------------------------------------------------------*/

PRIVATE void *
pop(LIST * list)
{
  void *data;

  data = lst_deletenext(list, LST_HEAD(list));
  return data;
}

/*---------------------------------------------------------------------------*/
/*auxiliary routines---------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

PRIVATE int
best_attainable_energy(STATE * state)
{
  /* evaluation of best possible energy attainable within remaining intervals */

  register int sum;
  INTERVAL *next;

  sum = state->partial_energy;  /* energy of already found elements */

  for (next = lst_first(state->Intervals); next; next = lst_next(next))
    {
      if (next->array_flag == 0)
        sum += (circular) ? Fc : f5[next->j];
      else if (next->array_flag == 1)
        sum += fML[indx[next->j] + next->i];
      else if (next->array_flag == 2)
        sum += c[indx[next->j] + next->i];
      else if (next->array_flag == 3)
        sum += fM1[indx[next->j] + next->i];
      else if (next->array_flag == 4)
        sum += fc[next->i];
      else if (next->array_flag == 5)
        sum += fc[next->j];
    }

  return sum;
}

/*---------------------------------------------------------------------------*/

PRIVATE void
push_back(STATE * state)
{
  push(Stack, copy_state(state));
  return;
}

/*---------------------------------------------------------------------------*/

PRIVATE char*
get_structure(STATE * state)
{
  char* structure;

  structure = strdup(state->structure);
  return structure;
}

/*---------------------------------------------------------------------------*/
PRIVATE int
compare(const void *solution1, const void *solution2)
{
  if (((SOLUTION *) solution1)->energy > ((SOLUTION *) solution2)->energy)
    return 1;
  if (((SOLUTION *) solution1)->energy < ((SOLUTION *) solution2)->energy)
    return -1;
  return strcmp(((SOLUTION *) solution1)->structure,
                ((SOLUTION *) solution2)->structure);
}

/*---------------------------------------------------------------------------*/

PRIVATE void make_output(SOLUTION *SL, FILE *fp)  /* prints stuff */
{
  SOLUTION *sol;

  for (sol = SL; sol->structure!=NULL; sol++)
    if (cut_point<0) fprintf(fp, "%s %6.2f\n", sol->structure, sol->energy);
    else {
      char *tStruc;
      tStruc=costring(sol->structure);
      fprintf(fp, "%s %6.2f\n", tStruc, sol->energy);
      free(tStruc);
    }
}

/*---------------------------------------------------------------------------*/
/* start of subopt backtracking ---------------------------------------------*/
/*---------------------------------------------------------------------------*/

PUBLIC SOLUTION *subopt(char *seq, char *structure, int delta, FILE *fp){
  return subopt_par(seq, structure, NULL, delta, fold_constrained, 0, fp);
}

PUBLIC SOLUTION *subopt_circ(char *seq, char *structure, int delta, FILE *fp){
  return subopt_par(seq, structure, NULL, delta, fold_constrained, 1, fp);
}

PUBLIC SOLUTION *subopt_par(char *seq,
                            char *structure,
                            paramT *parameters,
                            int delta,
                            int is_constrained,
                            int is_circular,
                            FILE *fp){

  STATE         *state;
  LIST          *Intervals;
  INTERVAL      *interval;
  SOLUTION      *SolutionList;
  unsigned long max_sol, n_sol;
  int           maxlevel, count, partial_energy, old_dangles, logML, dangle_model;
  double        structure_energy, min_en, eprint;
  char          *struc;

  max_sol             = 128;
  n_sol               = 0;
  sequence            = seq;
  length              = strlen(sequence);
  circular            = is_circular;
  struct_constrained  = is_constrained;

  struc = (char *) space(sizeof(char)*(length+1));
  if (struct_constrained) strncpy(struc, structure, length);

  /* do mfe folding to get fill arrays and get ground state energy  */
  /* in case dangles is neither 0 or 2, set dangles=2 while folding */

  if(P) free(P);
  if(parameters){
    P = get_parameter_copy(parameters);
  } else {
    model_detailsT md;
    set_model_details(&md);
    P = get_scaled_parameters(temperature, md);
  }

  logML       = P->model_details.logML;
  old_dangles = dangle_model = P->model_details.dangles;

  /* temporarily set dangles to 2 if necessary */
  if((P->model_details.dangles != 0) && (P->model_details.dangles != 2))
    P->model_details.dangles = 2;


  turn = (cut_point<0) ? 3 : 0;
  uniq_ML = 1;
  if(circular){
    min_en = fold_par(sequence, struc, P, struct_constrained, circular);
    export_circfold_arrays(&Fc, &FcH, &FcI, &FcM, &fM2, &f5, &c, &fML, &fM1, &indx, &ptype);
    /* restore dangle model */
    P->model_details.dangles = old_dangles;
    /* re-evaluate in case we're using logML etc */
    min_en = energy_of_circ_struct_par(sequence, struc, P, 0);
  } else {
    min_en = cofold_par(sequence, struc, P, struct_constrained);
    export_cofold_arrays(&f5, &c, &fML, &fM1, &fc, &indx, &ptype);
    /* restore dangle model */
    P->model_details.dangles = old_dangles;
    /* re-evaluate in case we're using logML etc */
    min_en = energy_of_struct_par(sequence, struc, P, 0);
  }

  free(struc);
  eprint = print_energy + min_en;
  if (fp) {
    char *SeQ;
    SeQ=costring(sequence);
    fprintf(fp, "%s %6d %6d\n", SeQ, (int) (-0.1+100*min_en), delta);
    free(SeQ);
  }
  make_pair_matrix();
  S   = encode_sequence(sequence, 0);
  S1  = encode_sequence(sequence, 1);

  /* Initialize ------------------------------------------------------------ */

  maxlevel = 0;
  count = 0;
  partial_energy = 0;

  /* Initialize the stack ------------------------------------------------- */

  minimal_energy = (circular) ? Fc : f5[length];
  threshold = minimal_energy + delta;
  if(threshold > INF){
    warn_user("energy range too high, limiting to reasonable value");
    threshold = INF-EMAX;
  }
  Stack = make_list();                                                   /* anchor */
  Intervals = make_list();                                   /* initial state: */
  interval = make_interval(1, length, 0);          /* interval [1,length,0] */
  push(Intervals, interval);
  state = make_state(Intervals, NULL, partial_energy,0);
  /* state->best_energy = minimal_energy; */
  push(Stack, state);

  /* SolutionList stores the suboptimal structures found */

  SolutionList = (SOLUTION *) space(max_sol*sizeof(SOLUTION));

  /* end initialize ------------------------------------------------------- */


  while (1) {                    /* forever, til nothing remains on stack */

    maxlevel = (Stack->count > maxlevel ? Stack->count : maxlevel);

    if (LST_EMPTY (Stack))                   /* we are done! clean up and quit */
      {
        /* fprintf(stderr, "maxlevel: %d\n", maxlevel); */

        lst_kill(Stack, free_state_node);

        SolutionList[n_sol].structure = NULL; /* NULL terminate list */

        if (subopt_sorted) {
          /* sort structures by energy */
          qsort(SolutionList, n_sol, sizeof(SOLUTION), compare);

          if (fp) make_output(SolutionList, fp);
        }

        break;
      }

    /* pop the last element ---------------------------------------------- */

    state = pop(Stack);                       /* current state to work with */

    if (LST_EMPTY(state->Intervals))
      {
        int e;
        /* state has no intervals left: we got a solution */

        count++;
        structure = get_structure(state);
        structure_energy = state->partial_energy / 100.;

#ifdef CHECK_ENERGY
        structure_energy = (circular) ? energy_of_circ_struct_par(sequence, structure, P, 0) : energy_of_struct_par(sequence, structure, P, 0);

        if (!logML)
          if ((double) (state->partial_energy / 100.) != structure_energy) {
            fprintf(stderr, "%s %6.2f %6.2f\n", structure,
                    state->partial_energy / 100., structure_energy );
            exit(1);
          }
#endif
        if (logML || (dangle_model==1) || (dangle_model==3)) { /* recalc energy */
          structure_energy = (circular) ? energy_of_circ_struct_par(sequence, structure, P, 0) : energy_of_struct_par(sequence, structure, P, 0);
        }

        e = (int) ((structure_energy-min_en)*10. + 0.1); /* avoid rounding errors */
        if (e>MAXDOS) e=MAXDOS;
        density_of_states[e]++;
        if (structure_energy>eprint) {
          free(structure);
        } else {
          if (!subopt_sorted && fp) {
            /* print and forget */
            if (cut_point<0)
              fprintf(fp, "%s %6.2f\n", structure, structure_energy);
            else {
              char * outstruc;
              /*make ampersand seperated output if 2 sequences*/
              outstruc=costring(structure);
              fprintf(fp, "%s %6.2f\n", outstruc, structure_energy);
              free(outstruc);
            }
            free(structure);
          }
          else {
            /* store solution */
            if (n_sol+1 == max_sol) {
              max_sol *= 2;
              SolutionList = (SOLUTION *)
                xrealloc(SolutionList, max_sol*sizeof(SOLUTION));
            }
            SolutionList[n_sol].energy =  structure_energy;
            SolutionList[n_sol++].structure = structure;
          }
        }
      }
    else {
      /* get (and remove) next interval of state to analyze */

      interval = pop(state->Intervals);

      scan_interval(interval->i, interval->j, interval->array_flag, state);

      free_interval_node(interval);        /* free the current interval */
    }

    free_state_node(state);                     /* free the current state */
  } /* end of while (1) */

  /* free arrays left over from cofold() */
  free(S); free(S1);
  (circular) ? free_arrays():free_co_arrays();
  if (fp) { /* we've printed everything -- free solutions */
    SOLUTION *sol;
    for (sol=SolutionList; sol->structure != NULL; sol++)
      free(sol->structure);
    free(SolutionList);
    SolutionList = NULL;
  }

  return SolutionList;
}


PRIVATE void
scan_interval(int i, int j, int array_flag, STATE * state)
{
  /* real backtrack routine */

  /* array_flag = 0:  trace back in f5-array  */
  /* array_flag = 1:  trace back in fML-array */
  /* array_flag = 2:  trace back in repeat()  */
  /* array_flag = 3:  trace back in fM1-array */

  STATE *new_state, *temp_state;
  INTERVAL *new_interval;
  register int k, fi, cij;
  register int type;
  register int dangle_model = P->model_details.dangles;
  register int noGUclosure  = P->model_details.noGUclosure;
  register int noLP         = P->model_details.noLP;

  best_energy = best_attainable_energy(state);  /* .. on remaining intervals */
  nopush = true;

  if ((i > 1) && (!array_flag))
    nrerror ("Error while backtracking!");

  if (j < i + turn + 1 && SAME_STRAND(i,j)) { /* minimal structure element */
    if (nopush)
      push_back(state);
    return;
  }

  /* 13131313131313131313131313131313131313131313131313131313131313131313131 */
  if (array_flag == 3 || array_flag == 1) {
    /* array_flag = 3: interval i,j was generated during */
    /*                 a multiloop decomposition using array fM1 in repeat() */
    /*                 or in this block */

    /* array_flag = 1: interval i,j was generated from a */
    /*                 stack, bulge, or internal loop in repeat() */
    /*                 or in this block */

    if (array_flag == 3)
      fi = fM1[indx[j-1] + i] + P->MLbase;
    else
      fi = fML[indx[j-1] + i] + P->MLbase;

    if ((fi + best_energy <= threshold)&&(SAME_STRAND(j-1,j))) {
      /* no basepair, nibbling of 3'-end */

      new_state = copy_state(state);
      new_interval = make_interval(i, j-1, array_flag);
      push(new_state->Intervals, new_interval);
      new_state->partial_energy += P->MLbase;
      /* new_state->best_energy = fi + best_energy; */
      push(Stack, new_state);
    }

    type = ptype[indx[j]+i];

    if (type) { /* i,j may pair */

      if(dangle_model)
        element_energy = E_MLstem(type,
                                  (((i > 1)&&(SAME_STRAND(i-1,i))) || circular)       ? S1[i-1] : -1,
                                  (((j < length)&&(SAME_STRAND(j,j+1))) || circular)  ? S1[j+1] : -1,
                                  P);
      else
        element_energy = E_MLstem(type, -1, -1, P);

      cij = c[indx[j] + i] + element_energy;
      if (cij + best_energy <= threshold)
        repeat(i, j, state, element_energy, 0);
    }
  }                                   /* array_flag == 3 || array_flag == 1 */

  /* 11111111111111111111111111111111111111111111111111111111111111111111111 */

  if (array_flag == 1) {
    /* array_flag = 1:                   interval i,j was generated from a */
    /*                          stack, bulge, or internal loop in repeat() */
    /*                          or in this block */

    int stopp;
    if ((SAME_STRAND(i-1,i))&&(SAME_STRAND(j,j+1))) { /*backtrack in FML only if multiloop is possible*/
      for ( k = i+turn+1 ; k <= j-1-turn ; k++) {
        /* Multiloop decomposition if i,j contains more than 1 stack */

        type = ptype[indx[j]+k+1];
        if (type==0) continue;

        if(dangle_model)
          element_energy = E_MLstem(type,
                                    (SAME_STRAND(i-1,i)) ? S1[k] : -1,
                                    (SAME_STRAND(j,j+1)) ? S1[j+1] : -1,
                                    P);
        else
          element_energy = E_MLstem(type, -1, -1, P);

        if (SAME_STRAND(k,k+1)) {
          if (fML[indx[k]+i] + c[indx[j] + k+1] +
              element_energy + best_energy <= threshold) {
            temp_state = copy_state (state);
            new_interval = make_interval (i, k, 1);
            push (temp_state->Intervals, new_interval);
            repeat(k+1, j, temp_state, element_energy, fML[indx[k]+i]);
            free_state_node(temp_state);
          }
        }
      }
    }

    stopp=(cut_point>0)? (cut_point-2):(length); /*if cut_point -1: k on cut, => no ml*/
    stopp=MIN2(stopp, j-1-turn);
    if (i>cut_point) stopp=j-1-turn;
    else if (i==cut_point) stopp=0;   /*not a multi loop*/
    for (k = i ; k <= stopp; k++) {
      /* Multiloop decomposition if i,j contains only 1 stack */
      type = ptype[indx[j]+k+1];
      if (type==0) continue;

      if(dangle_model)
        element_energy = E_MLstem(type,
                                  (SAME_STRAND(k-1,k)) ? S1[k] : -1,
                                  (SAME_STRAND(j,j+1)) ? S1[j+1] : -1,
                                  P);
      else
        element_energy = E_MLstem(type, -1, -1, P);

      element_energy += P->MLbase*(k-i+1);

      if (c[indx[j]+k+1] + element_energy + best_energy <= threshold)
        repeat(k+1, j, state, element_energy, 0);
    }
  }                                                    /* array_flag == 1 */

  /* 2222222222222222222222222222222222222222222222222222222222222222222222 */

  if (array_flag == 2)
    {
      /* array_flag = 2:                  interval i,j was generated from a */
      /* stack, bulge, or internal loop in repeat() */

      repeat(i, j, state, 0, 0);

      if (nopush){
        if (!noLP){
          fprintf(stderr, "%d,%d", i, j);
          fprintf(stderr, "Oops, no solution in repeat!\n");
        }
      }
      return;
    }

  /* 0000000000000000000000000000000000000000000000000000000000000000000000 */

  if ((array_flag == 0) && !circular)
    {
      /* array_flag = 0:                        interval i,j was found while */
      /* tracing back through f5-array and c-array */
      /* or within this block */

      if (f5[j-1] + best_energy <= threshold) {
        /* no basepair, nibbling of 3'-end */

        new_state = copy_state(state);
        new_interval = make_interval(i, j-1 , 0);
        push(new_state->Intervals, new_interval);
        /* new_state->best_energy = f5[j-1] + best_energy; */
        push(Stack, new_state);
      }

      for (k = j-turn-1; k > 1; k--) {

        type = ptype[indx[j]+k];
        if (type==0)   continue;

        /* k and j pair */
        if(dangle_model)
          element_energy = E_ExtLoop(type,
                                    (SAME_STRAND(k-1,k)) ? S1[k-1] : -1,
                                    ((j < length)&&(SAME_STRAND(j,j+1))) ? S1[j+1] : -1,
                                    P);
        else
          element_energy = E_ExtLoop(type, -1, -1, P);

        if (!(SAME_STRAND(k,j)))/*&&(state->is_duplex==0))*/ {
          element_energy+=P->DuplexInit;
          /*state->is_duplex=1;*/
        }

        if (f5[k-1] + c[indx[j]+k] + element_energy + best_energy <= threshold)
          {
            temp_state = copy_state(state);
            new_interval = make_interval(1, k-1, 0);
            push(temp_state->Intervals, new_interval);
            repeat(k, j, temp_state, element_energy, f5[k-1]);
            free_state_node(temp_state);
          }
      }
      type = ptype[indx[j]+1];
      if (type) {
        if (dangle_model && (j < length)&&(SAME_STRAND(j,j+1)))
          element_energy = E_ExtLoop(type, -1, S1[j+1], P);
        else
          element_energy = E_ExtLoop(type, -1, -1, P);

        if (!(SAME_STRAND(1,j))) element_energy+=P->DuplexInit;

        if (c[indx[j]+1] + element_energy + best_energy <= threshold)
          repeat(1, j, state, element_energy, 0);
      }
    }/* end array_flag == 0 && !circular*/
  /* or do we subopt circular? */
  else if(array_flag == 0){
    int k, l, p, q;
    /* if we've done everything right, we will never reach this case more than once         */
    /* right after the initilization of the stack with ([1,n], empty, 0)                                                 */
    /* lets check, if we can have an open chain without breaking the threshold        */
    /* this is an ugly work-arround cause in case of an open chain we do not have to  */
    /* backtrack anything further...                                                  */
    if(0 <= threshold){
      new_state = copy_state(state);
      new_interval = make_interval(1,2,0);
      push(new_state->Intervals, new_interval);
      new_state->partial_energy = 0;
      push(Stack, new_state);
    }
    /* ok, lets check if we can do an exterior hairpin without breaking the threshold */
    /* best energy should be 0 if we are here                                                                                                                                                                 */
    if(FcH + best_energy <= threshold){
      /* lets search for all exterior hairpin cases, that fit into our threshold barrier */
      /* we use index k,l to avoid confusion with i,j index of our state...                                                         */
      /* if we reach here, i should be 1 and j should be n respectively                                                                         */
      for(k=i; k<j; k++)
        for (l=k+turn+1; l <= j; l++){
          int kl, type, u, new_c, tmpE, no_close;
          u = j-l + k-1;        /* get the hairpin loop length */
          if(u<turn) continue;

          kl = indx[l]+k;        /* just confusing these indices ;-) */
          type = ptype[kl];
          no_close = ((type==3)||(type==4))&&noGUclosure;
          type=rtype[type];
          if (!type) continue;
          if (no_close) new_c = FORBIDDEN;
          else{
            /* now lets have a look at the hairpin energy */
            char loopseq[10];
            if (u<7){
              strcpy(loopseq , sequence+l-1);
              strncat(loopseq, sequence, k);
            }
            tmpE = E_Hairpin(u, type, S1[l+1], S1[k-1], loopseq, P);
          }
          if(c[kl] + tmpE + best_energy <= threshold){
            /* what we really have to do is something like this, isn't it? */
            /* we have to create a new state, with interval [k,l], then we */
            /* add our loop energy as initial energy of this state and put */
            /* the state onto the stack R... for further refinement...         */
            /* we also denote this new interval to be scanned in C         */
            new_state = copy_state(state);
            new_interval = make_interval(k,l,2);
            push(new_state->Intervals, new_interval);
            /* hopefully we add this energy in the right way... */
            new_state->partial_energy += tmpE;
            push(Stack, new_state);
          }
        }
    }

    /* now lets see, if we can do an exterior interior loop without breaking the threshold */
    if(FcI + best_energy <= threshold){
      /* now we search for our exterior interior loop possibilities */
      for(k=i; k<j; k++)
        for (l=k+turn+1; l <= j; l++){
          int kl, type, tmpE, no_close;

          kl = indx[l]+k;        /* just confusing these indices ;-) */
          type = ptype[kl];
          no_close = ((type==3)||(type==4))&&noGUclosure;
          type=rtype[type];
          if (!type) continue;

          for (p = l+1; p < j ; p++){
            int u1, qmin;
            u1 = p-l-1;
            if (u1+k-1>MAXLOOP) break;
            qmin = u1+k-1+j-MAXLOOP;
            if(qmin<p+turn+1) qmin = p+turn+1;
            for(q = qmin; q <=j; q++){
              int u2, type_2;
              type_2 = rtype[ptype[indx[q]+p]];
              if(!type_2) continue;
              u2 = k-1 + j-q;
              if(u1+u2>MAXLOOP) continue;
              tmpE = E_IntLoop(u1, u2, type, type_2, S1[l+1], S1[k-1], S1[p-1], S1[q+1], P);
              if(c[kl] + c[indx[q]+p] + tmpE + best_energy <= threshold){
                /* ok, similar to the hairpin stuff, we add new states onto the stack R */
                /* but in contrast to the hairpin decomposition, we have to add two new */
                /* intervals, enclosed by k,l and p,q respectively and we also have to */
                /* add the partial energy, that comes from the exterior interior loop   */
                new_state = copy_state(state);
                new_interval = make_interval(k, l, 2);
                push(new_state->Intervals, new_interval);
                new_interval = make_interval(p,q,2);
                push(new_state->Intervals, new_interval);
                new_state->partial_energy += tmpE;
                push(Stack, new_state);
              }
            }
          }
        }
    }

    /* and last but not least, we have a look, if we can do an exterior multiloop within the energy threshold */
    if(FcM <= threshold){
      /* this decomposition will be somehow more complicated...so lets see what we do here...           */
      /* first we want to find out which split inidices we can use without exceeding the threshold */
      int tmpE2;
      for (k=turn+1; k<j-2*turn; k++){
        tmpE2 = fML[indx[k]+1]+fM2[k+1]+P->MLclosing;
        if(tmpE2 + best_energy <= threshold){
          /* grmpfh, we have found a possible split index k so we have to split fM2 and fML now */
          /* lets do it first in fM2 anyway */
          for(l=k+turn+2; l<j-turn-1; l++){
            tmpE2 = fM1[indx[l]+k+1] + fM1[indx[j]+l+1];
            if(tmpE2 + fML[indx[k]+1] + P->MLclosing <= threshold){
              /* we've (hopefully) found a valid decomposition of fM2 and therefor we have all */
              /* three intervals for our new state to be pushed on stack R */
              new_state = copy_state(state);

              /* first interval leads for search in fML array */
              new_interval = make_interval(1, k, 1);
              push(new_state->Intervals, new_interval);

              /* next, we have the first interval that has to be traced in fM1 */
              new_interval = make_interval(k+1, l, 3);
              push(new_state->Intervals, new_interval);

              /* and the last of our three intervals is also one to be traced within fM1 array... */
              new_interval = make_interval(l+1, j, 3);
              push(new_state->Intervals, new_interval);

              /* mmh, we add the energy for closing the multiloop now... */
              new_state->partial_energy += P->MLclosing;
              /* next we push our state onto the R stack */
              push(Stack, new_state);

            }
            /* else we search further... */
          }

          /* ok, we have to decompose fML now... */
        }
      }
    }
  }        /* thats all folks for the circular case... */

  /* 44444444444444444444444444444444444444444444444444444444444444 */
  if (array_flag == 4) {
    /* array_flag = 4:                        interval i,j was found while */
    /* tracing back through fc-array smaller than than cut_point*/
    /* or within this block */

    if (fc[i+1] + best_energy <= threshold) {
      /* no basepair, nibbling of 5'-end */
      new_state = copy_state(state);
      new_interval = make_interval(i+1, j , 4);
      push(new_state->Intervals, new_interval);
      push(Stack, new_state);
    }

    for (k = i+TURN+1; k < j; k++) {

      type = ptype[indx[k]+i];
      if (type==0)   continue;

      /* k and j pair */
      if (dangle_model)
        element_energy = E_ExtLoop(type, (i > 1) ? S1[i-1]: -1, S1[k+1], P);
      else  /* no dangles */
        element_energy = E_ExtLoop(type, -1, -1, P);

      if (fc[k+1] + c[indx[k]+i] + element_energy + best_energy <= threshold) {
        temp_state = copy_state(state);
        new_interval = make_interval(k+1,j, 4);
        push(temp_state->Intervals, new_interval);
        repeat(i, k, temp_state, element_energy, fc[k+1]);
        free_state_node(temp_state);
      }
    }
    type = ptype[indx[j]+i];
    if (type) {
      if (dangle_model)
        element_energy = E_ExtLoop(type, (i>1) ? S1[i-1] : -1, -1, P);
      else
        element_energy = E_ExtLoop(type, -1, -1, P);

      if (c[indx[cut_point-1]+i] + element_energy + best_energy <= threshold)
        repeat(i, cut_point-1, state, element_energy, 0);
    }
  } /* array_flag == 4 */

  /*55555555555555555555555555555555555555555555555555555555555555555555555*/
  if (array_flag == 5) {
    /* array_flag = 5:  interval cut_point=i,j was found while  */
    /* tracing back through fc-array greater than cut_point     */
    /* or within this block                                     */

    if (fc[j-1] + best_energy <= threshold) {
      /* no basepair, nibbling of 3'-end */
      new_state = copy_state(state);
      new_interval = make_interval(i, j-1 , 5);
      push(new_state->Intervals, new_interval);
      push(Stack, new_state);
    }

    for (k = j-TURN-1; k > i; k--) {

      type = ptype[indx[j]+k];
      if (type==0)   continue;
      element_energy = 0;

      /* k and j pair */
      if (dangle_model)
        element_energy = E_ExtLoop(type, S1[k-1], (j < length) ? S1[j+1] : -1, P);
      else
        element_energy = E_ExtLoop(type, -1, -1, P);

      if (fc[k-1] + c[indx[j]+k] + element_energy + best_energy <= threshold) {
        temp_state = copy_state(state);
        new_interval = make_interval(i, k-1, 5);
        push(temp_state->Intervals, new_interval);
        repeat(k, j, temp_state, element_energy, fc[k-1]);
        free_state_node(temp_state);
      }
    }
    type = ptype[indx[j]+i];
    if (type) {
      if(dangle_model)
        element_energy = E_ExtLoop(type, -1, (j<length) ? S1[j+1] : -1, P);

      if (c[indx[j]+cut_point] + element_energy + best_energy <= threshold)
        repeat(cut_point, j, state, element_energy, 0);
    }
  } /* array_flag == 5 */
  if (nopush)
    push_back(state);
  return;
}

/*---------------------------------------------------------------------------*/

PRIVATE void
repeat(int i, int j, STATE * state, int part_energy, int temp_energy)
{
  /* routine to find stacks, bulges, internal loops and  multiloops */
  /* within interval closed by basepair i,j */

  STATE *new_state;
  INTERVAL *new_interval;

  register int  k, p, q, energy, new;
  register int  mm;
  register int  no_close, no_close_2, type, type_2;
  int           rt;
  int           dangle_model  = P->model_details.dangles;
  int           noLP          = P->model_details.noLP;
  int           noGUclosure   = P->model_details.noGUclosure;

  no_close_2 = 0;

  type = ptype[indx[j]+i];
  if (type==0) fprintf(stderr, "repeat: Warning: %d %d can't pair\n", i,j);

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

  if (noLP) /* always consider the structure with additional stack */
    if ((i+turn+2<j) && ((type_2 = ptype[indx[j-1]+i+1]))) {
      new_state = copy_state(state);
      make_pair(i, j, new_state);
      make_pair(i+1, j-1, new_state);
      new_interval = make_interval(i+1, j-1, 2);
      push(new_state->Intervals, new_interval);
      if(SAME_STRAND(i,i+1) && SAME_STRAND(j-1,j))
        energy = E_IntLoop(0, 0, type, rtype[type_2],S1[i+1],S1[j-1],S1[i+1],S1[j-1], P);

      new_state->partial_energy += part_energy;
      new_state->partial_energy += energy;
      /* new_state->best_energy = new + best_energy; */
      push(Stack, new_state);
      if (i==1 || state->structure[i-2]!='('  || state->structure[j]!=')')
        /* adding a stack is the only possible structure */
        return;
    }

  best_energy += part_energy; /* energy of current structural element */
  best_energy += temp_energy; /* energy from unpushed interval */

  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 = j - 1; q >= minq; q--) {
      if ((noLP) && (p==i+1) && (q==j-1)) continue;

      type_2 = ptype[indx[q]+p];
      if (type_2==0) continue;

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

      if (SAME_STRAND(i,p) && SAME_STRAND(q,j)) {
        energy = E_IntLoop(p-i-1, j-q-1, type, rtype[type_2],
                            S1[i+1],S1[j-1],S1[p-1],S1[q+1], P);

        new = energy + c[indx[q]+p];

        if (new + best_energy <= threshold) {
          /* stack, bulge, or interior loop */

          new_state = copy_state(state);
          make_pair(i, j, new_state);
          make_pair(p, q, new_state);

          new_interval = make_interval(p, q, 2);
          push(new_state->Intervals, new_interval);
          new_state->partial_energy += part_energy;
          new_state->partial_energy += energy;
          /* new_state->best_energy = new + best_energy; */
          push(Stack, new_state);
        }
      }/*end of if block */
    } /* end of q-loop */
  } /* end of p-loop */

  if (!SAME_STRAND(i,j)) { /*look in fc*/
    rt = rtype[type];
    element_energy=0;
    if (dangle_model)
      element_energy = E_ExtLoop(rt, (SAME_STRAND(j-1,j)) ? S1[j-1] : -1, (SAME_STRAND(i,i+1)) ? S1[i+1] : -1, P);
    else
      element_energy = E_ExtLoop(rt, -1, -1, P);

    if (fc[i+1] + fc[j-1] +element_energy + best_energy  <= threshold)
      {
        INTERVAL *interval1, *interval2;

        new_state = copy_state(state);
        interval1 = make_interval(i+1, cut_point-1, 4);
        interval2 = make_interval(cut_point, j-1, 5);
        if (cut_point-i < j-cut_point) { /* push larger interval first */
          push(new_state->Intervals, interval1);
          push(new_state->Intervals, interval2);
        } else {
          push(new_state->Intervals, interval2);
          push(new_state->Intervals, interval1);
        }
        make_pair(i, j, new_state);
        new_state->partial_energy += part_energy;
        new_state->partial_energy += element_energy;
        push(Stack, new_state);
      }
  }

  mm = P->MLclosing;
  rt = rtype[type];

  for (k = i + 1 + turn; k <= j - 2 - turn; k++)  {
    /* multiloop decomposition */

    element_energy = mm;
    if (dangle_model)
      element_energy = E_MLstem(rt, S1[j-1], S1[i+1], P) + mm;
    else
      element_energy = E_MLstem(rt, -1, -1, P) + mm;

    if ((fML[indx[k] + i+1] + fM1[indx[j-1] + k+1] +
        element_energy + best_energy)  <= threshold)
      {
        INTERVAL *interval1, *interval2;

        new_state = copy_state(state);
        interval1 = make_interval(i+1, k, 1);
        interval2 = make_interval(k+1, j-1, 3);
        if (k-i+1 < j-k-2) { /* push larger interval first */
          push(new_state->Intervals, interval1);
          push(new_state->Intervals, interval2);
        } else {
          push(new_state->Intervals, interval2);
          push(new_state->Intervals, interval1);
        }
        make_pair(i, j, new_state);
        new_state->partial_energy += part_energy;
        new_state->partial_energy += element_energy;
        /* new_state->best_energy = fML[indx[k] + i+1] + fM1[indx[j-1] + k+1]
           + element_energy + best_energy; */
        push(Stack, new_state);
      }
  } /* end of k-loop */


  if (SAME_STRAND(i,j)) {
    if (no_close) element_energy = FORBIDDEN;
    else
      element_energy = E_Hairpin(j-i-1, type, S1[i+1], S1[j-1], sequence+i-1, P);
    if (element_energy + best_energy <= threshold) {
      /* hairpin structure */

      new_state = copy_state(state);
      make_pair(i, j, new_state);
      new_state->partial_energy += part_energy;
      new_state->partial_energy += element_energy;
      /* new_state->best_energy =
         hairpin[unpaired] + element_energy + best_energy; */
      push(Stack, new_state);
    }
  }

  best_energy -= part_energy;
  best_energy -= temp_energy;
  return;
}

PRIVATE char *costring(char *string)
{
  char *ctmp;
  int len;
  len = strlen(string);
  ctmp = (char *)space((len+2) * sizeof(char));
  /* first sequence */
  if (cut_point<=0) {
    (void) strncpy(ctmp, string, len);
    return ctmp;
  }
  (void) strncpy(ctmp, string, cut_point-1);
  /* spacer */
  ctmp[cut_point-1] = '&';
  /* second sequence */
  (void) strcat(ctmp, string+cut_point-1);

  return ctmp;
}

/*---------------------------------------------------------------------------*/
/* Well, that is the end!----------------------------------------------------*/
/*---------------------------------------------------------------------------*/