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

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

#include "ViennaRNA/pair_mat.h"
#include "ViennaRNA/fold_vars.h"
#include "ViennaRNA/fold.h"
#include "ViennaRNA/utils.h"

#include "ViennaRNA/move_set.h"


/* maximum degeneracy value - if degeneracy is greater than this, program segfaults */
#define MAX_DEGEN 100
#define MINGAP 3


#define bool int
#define true 1
#define false 0

/*
#################################
# PRIVATE DATA STRUCTURES       #
#################################
*/

/* internal struct with moves, sequence, degeneracy and options*/
typedef struct _Encoded {
  /* sequence*/
  short *s0;
  short *s1;

  const char  *seq;

  /* moves*/
  int   bp_left;
  int   bp_right;
  int   bp_left2;   /* if noLP is enabled (and for shift moves)*/
  int   bp_right2;

  /* options*/
  int noLP;
  int verbose_lvl;
  int first;
  int shift;

  /* degeneracy*/
  int begin_unpr;
  int begin_pr;
  int end_unpr;
  int end_pr;
  short *processed[MAX_DEGEN];
  short *unprocessed[MAX_DEGEN];
  int current_en;

  /* moves in random (needs to be freed afterwards)*/
  int *moves_from;
  int *moves_to;
  int num_moves;

  /* function for flooding */
  int (*funct) (struct_en*, struct_en*);


} Encoded;

/*
#################################
# PRIVATE VARIABLES             #
#################################
*/
PRIVATE int cnt_move = 0;

/*
#################################
# PRIVATE FUNCTION DECLARATIONS #
#################################
*/
PRIVATE int     compare(short *lhs, short *rhs);
PRIVATE int     find_min(short *arr[MAX_DEGEN], int begin, int end);
PRIVATE int     equals(const short *first, const short *second);
PRIVATE int     count_move(void);
PRIVATE int     lone_base(short *pt, int i);
PRIVATE void    free_degen(Encoded *Enc);
PRIVATE inline void do_move(short *pt, int bp_left, int bp_right);
PRIVATE int     update_deepest(Encoded *Enc, struct_en *str, struct_en *min);
PRIVATE int     deletions(Encoded *Enc, struct_en *str, struct_en *minim);
PRIVATE inline  bool compat(char a, char b);
PRIVATE inline  bool try_insert(const short *pt, const char *seq, int i, int j);
PRIVATE inline  bool try_insert_seq(const char *seq, int i, int j);
PRIVATE int     insertions(Encoded *Enc, struct_en *str, struct_en *minim);
PRIVATE int     shifts(Encoded *Enc, struct_en *str, struct_en *minim);
PRIVATE int     move_set(Encoded *Enc, struct_en *str);
PRIVATE void    construct_moves(Encoded *Enc, short *structure);
PRIVATE int     move_rset(Encoded *Enc, struct_en *str);
PRIVATE int     find_lone_pair(short* str);


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

PRIVATE int
compare(short *lhs, short *rhs){

  /* printf("%d ", (int)lhs[0]); */

  int i=1;
  char l=0,r=0;
  while (i<=lhs[0]) {
    l = (lhs[i]==0?'.':(lhs[i]<lhs[lhs[i]]?')':'('));
    r = (rhs[i]==0?'.':(rhs[i]<rhs[rhs[i]]?')':'('));
    if (l != r) break;
    i++;
  }

  return (i<=lhs[0] && l>r);
}

PRIVATE int
find_min(short *arr[MAX_DEGEN], int begin, int end){

  short *min = arr[begin];
  short min_num = begin;
  int i;

  for (i=begin+1; i<end; i++) {
    if (compare(arr[i], min)) {
      min = arr[i];
      min_num = i;
    }
  }
  return min_num;
}

PRIVATE int
equals(const short *first, const short *second){

  int i=1;
  while (i<=first[0] && first[i]==second[i]) {
    i++;
  }
  if (i>first[0]) return 1;
  else return 0;
}

PUBLIC void
copy_arr(short *dest, short *src){
  if (!src || !dest) {
    vrna_message_warning("Empty pointer in copying");
    return;
  }
  memcpy(dest, src, sizeof(short)*(src[0]+1));
}

PUBLIC short *
allocopy(short *src){
  short *res = (short*) vrna_alloc(sizeof(short)*(src[0]+1));
  copy_arr(res, src);
  return res;
}

PRIVATE int
count_move(void){

  return cnt_move;
}



/* frees all things allocated by degeneracy...*/
PRIVATE void
free_degen(Encoded *Enc){

  int i;
  for (i=Enc->begin_unpr; i<Enc->end_unpr; i++) {
    if (Enc->unprocessed[i]) {
      free(Enc->unprocessed[i]);
      Enc->unprocessed[i]=NULL;
    }
  }
  for (i=Enc->begin_pr; i<Enc->end_pr; i++) {
    if (Enc->processed[i]) {
      free(Enc->processed[i]);
      Enc->processed[i]=NULL;
    }
  }
  Enc->begin_pr=0;
  Enc->begin_unpr=0;
  Enc->end_pr=0;
  Enc->end_unpr=0;
}

PRIVATE inline void
do_move(short *pt, int bp_left, int bp_right){

  /* delete*/
  if (bp_left<0) {
    pt[-bp_left]=0;
    pt[-bp_right]=0;
  } else { /* insert*/
    pt[bp_left]=bp_right;
    pt[bp_right]=bp_left;
  }
}

/* done with all structures along the way to deepest*/
PRIVATE int
update_deepest(Encoded *Enc, struct_en *str, struct_en *min){

  /* apply move + get its energy*/
  int tmp_en;
  tmp_en = str->energy + energy_of_move_pt(str->structure, Enc->s0, Enc->s1, Enc->bp_left, Enc->bp_right);
  do_move(str->structure, Enc->bp_left, Enc->bp_right);
  if (Enc->bp_left2 != 0) {
    tmp_en += energy_of_move_pt(str->structure, Enc->s0, Enc->s1, Enc->bp_left2, Enc->bp_right2);
    do_move(str->structure, Enc->bp_left2, Enc->bp_right2);
  }
  int last_en = str->energy;
  str->energy = tmp_en;


  /* use f_point if we have it */
  if (Enc->funct) {
    int end = Enc->funct(str, min);

    /*  undo moves */
    if (Enc->bp_left2!=0) do_move(str->structure, -Enc->bp_left2, -Enc->bp_right2);
    do_move(str->structure, -Enc->bp_left, -Enc->bp_right);
    str->energy = last_en;
    Enc->bp_left=0;
    Enc->bp_right=0;
    Enc->bp_left2=0;
    Enc->bp_right2=0;

    return (end?1:0);
  }

  if (Enc->verbose_lvl>1) { fprintf(stderr, "  "); print_str(stderr, str->structure); fprintf(stderr, " %d\n", tmp_en); }

  /* better deepest*/
  if (tmp_en < min->energy) {
    min->energy = tmp_en;
    copy_arr(min->structure, str->structure);

    /* delete degeneracy*/
    free_degen(Enc);

    /* undo moves*/
    if (Enc->bp_left2!=0) do_move(str->structure, -Enc->bp_left2, -Enc->bp_right2);
    do_move(str->structure, -Enc->bp_left, -Enc->bp_right);
    str->energy = last_en;
    Enc->bp_left=0;
    Enc->bp_right=0;
    Enc->bp_left2=0;
    Enc->bp_right2=0;
    return 1;
  }

  /* degeneracy*/
  if ((str->energy == min->energy) && (Enc->current_en == min->energy)) {
    int found = 0;
    int i;
    for (i=Enc->begin_pr; i<Enc->end_pr; i++) {
      if (equals(Enc->processed[i], str->structure)) {
        found = 1;
        break;
      }
    }
    for (i=Enc->begin_unpr; !found && i<Enc->end_unpr; i++) {
      if (equals(Enc->unprocessed[i], str->structure)) {
        found = 1;
        break;
      }
    }

    if (!found) {
      /* print_stren(stderr, str); // fprintf(stderr, " %6.2f\n", str->energy); */
      Enc->unprocessed[Enc->end_unpr]=allocopy(str->structure);
      Enc->end_unpr++;
    }
  }

  /* undo moves*/
  if (Enc->bp_left2!=0) do_move(str->structure, -Enc->bp_left2, -Enc->bp_right2);
  do_move(str->structure, -Enc->bp_left, -Enc->bp_right);
  str->energy = last_en;
  Enc->bp_left=0;
  Enc->bp_right=0;
  Enc->bp_left2=0;
  Enc->bp_right2=0;
  return 0;
}


/* deletions move set*/
PRIVATE int
deletions(Encoded *Enc, struct_en *str, struct_en *minim){

  int cnt = 0;
  short *pt = str->structure;
  int len = pt[0];
  int i;

  for (i=1; i<=len; i++) {
    if (pt[i]>pt[pt[i]]) {  /* '('*/
      Enc->bp_left=-i;
      Enc->bp_right=-pt[i];

      /*if nolp enabled, make (maybe) 2nd delete*/
      if (Enc->noLP) {
        int lone = -1;
        if (lone_base(pt, i-1)) lone=i-1;
        else if (lone_base(pt, i+1)) lone=i+1;
        else if (lone_base(pt, pt[i]-1)) lone=pt[i]-1;
        else if (lone_base(pt, pt[i]+1)) lone=pt[i]+1;

        /* check*/
        if (lone != -1 && (pt[lone]==0 || pt[pt[lone]]==0)) {
          vrna_message_warning("pt[%d(or %d)]!=\'.\'", lone, pt[lone]);
        }

        if (lone != -1) {
          Enc->bp_left2=-lone-1;
          Enc->bp_right2=-pt[lone]-1;
        }
        if (!lone_base(pt, pt[lone]-1) && !lone_base(pt, pt[lone]+1)) {
          cnt += update_deepest(Enc, str, minim);
          /* in case useFirst is on and structure is found, end*/
          if (Enc->first && cnt > 0) return cnt;
        }
      } else {  /* nolp not enabled*/
        cnt += update_deepest(Enc, str, minim);
        /* in case useFirst is on and structure is found, end*/
        if (Enc->first && cnt > 0) return cnt;
      }
    }
  }
  return cnt;
}

  /* compatible base pair?*/
PRIVATE inline bool
compat(char a, char b){

  if (a=='A' && b=='U') return true;
  if (a=='C' && b=='G') return true;
  if (a=='G' && b=='U') return true;
  if (a=='U' && b=='A') return true;
  if (a=='G' && b=='C') return true;
  if (a=='U' && b=='G') return true;
  /* and with T's*/
  if (a=='A' && b=='T') return true;
  if (a=='T' && b=='A') return true;
  if (a=='G' && b=='T') return true;
  if (a=='T' && b=='G') return true;
  return false;
}

/* try insert base pair (i,j)*/
PRIVATE inline bool
try_insert(const short *pt, const char *seq, int i, int j){

  if (i<=0 || j<=0 || i>pt[0] || j>pt[0]) return false;
  return (j-i>MINGAP && pt[j]==0 && pt[i]==0 && compat(seq[i-1], seq[j-1]));
}

/*  try insert base pair (i,j) */
PRIVATE inline bool
try_insert_seq(const char *seq, int i, int j){
  if (i<=0 || j<=0) return false;
  return (j-i>MINGAP && compat(seq[i-1], seq[j-1]));
}

/* insertions move set*/
PRIVATE int
insertions(Encoded *Enc, struct_en *str, struct_en *minim){

  int cnt = 0;
  short *pt = str->structure;
  int len = pt[0];
  int i,j;

  for (i=1; i<=len; i++) {
    if (pt[i]==0) {
      for (j=i+1; j<=len; j++) {
        /* end if found closing bracket*/
        if (pt[j]!=0 && pt[j]<j) break;  /*')'*/
        if (pt[j]!=0 && pt[j]>j) {       /*'('*/
          j = pt[j];
          continue;
        }
        /* if conditions are met, do insert*/
        if (try_insert(pt, Enc->seq, i, j)) {
          Enc->bp_left=i;
          Enc->bp_right=j;

          if (Enc->noLP) {
            /* if lone bases occur, try inserting one another base*/
            if (lone_base(pt, i) || lone_base(pt, j)) {
              /* inside*/
              if (try_insert(pt, Enc->seq, i+1, j-1)) {
                Enc->bp_left2=i+1;
                Enc->bp_right2=j-1;
                cnt += update_deepest(Enc, str, minim);
                /* in case useFirst is on and structure is found, end*/
                if (Enc->first && cnt > 0) return cnt;
              } else  /*outside*/
              if (try_insert(pt, Enc->seq, i-1, j+1)) {
                Enc->bp_left2=i-1;
                Enc->bp_right2=j+1;
                cnt += update_deepest(Enc, str, minim);
                /* in case useFirst is on and structure is found, end*/
                if (Enc->first && cnt > 0) return cnt;
              }
            } else {
              cnt += update_deepest(Enc, str, minim);
              /* in case useFirst is on and structure is found, end*/
              if (Enc->first && cnt > 0) return cnt;
            }
          } else {
            cnt += update_deepest(Enc, str, minim);
            /* in case useFirst is on and structure is found, end*/
            if (Enc->first && cnt > 0) return cnt;
          }
        }
      }
    }
  }
  return cnt;
}

/*shift move set*/
PRIVATE int
shifts(Encoded *Enc, struct_en *str, struct_en *minim){

  int cnt = 0;
  int brack_num = 0;
  short *pt = str->structure;
  int len = pt[0];
  int i, k;

  for (i=1; i<=len; i++) {
    if (pt[i]!=0 && pt[i]>i) {  /*'('*/
      int j=pt[i];

      /* outer switch left*/
      if (Enc->verbose_lvl>1)
        vrna_message_info(stderr, "%2d bracket %2d position, outer switch left", brack_num+1, i);
      for (k=i-1; k>0; k--) {
        if (pt[k]!=0 && pt[k]>k/*'('*/) break;
        if (pt[k]!=0 && pt[k]<k/*')'*/) {
          k = pt[k];
          continue;
        }
        /* checks*/
        if (pt[k]!=0) {
          vrna_message_warning("\'%c\'should be \'.\' at pos %d!", pt[k], k);
        }

        /* switch (i,j) to (k,j)*/
        if (j-k>MINGAP && compat(Enc->seq[k-1], Enc->seq[j-1])) {
          Enc->bp_left=-i;
          Enc->bp_right=-j;
          Enc->bp_left2=k;
          Enc->bp_right2=j;
          cnt += update_deepest(Enc, str, minim);
          /* in case useFirst is on and structure is found, end*/
          if (Enc->first && cnt > 0) return cnt;
        }

        /* switch (i,j) to (k,i)*/
        if (i-k>MINGAP && compat(Enc->seq[i-1], Enc->seq[k-1])) {
          Enc->bp_left=-i;
          Enc->bp_right=-j;
          Enc->bp_left2=k;
          Enc->bp_right2=i;
          cnt += update_deepest(Enc, str, minim);
          /* in case useFirst is on and structure is found, end*/
          if (Enc->first && cnt > 0) return cnt;

        }
      }

      /* outer switch right*/
      if (Enc->verbose_lvl>1)
        vrna_message_info(stderr, "%2d bracket %2d position, outer switch right", brack_num+1, i);
      for (k=j+1; k<=len; k++) {
        if (pt[k]!=0 && pt[k]<k/*')'*/) break;
        if (pt[k]!=0 && pt[k]>k/*'('*/) {
          k = pt[k];
          continue;
        }

        /* check*/
        if (pt[k]!=0) {
          vrna_message_warning("\'%c\'should be \'.\' at pos %d!", pt[k], k);
        }
        /* switch (i,j) to (i,k)*/
        if (k-i>MINGAP && compat(Enc->seq[i-1], Enc->seq[k-1])) {
          Enc->bp_left=-i;
          Enc->bp_right=-j;
          Enc->bp_left2=i;
          Enc->bp_right2=k;
          cnt += update_deepest(Enc, str, minim);
          /* in case useFirst is on and structure is found, end*/
          if (Enc->first && cnt > 0) return cnt;
        }
        /* switch (i,j) to (j,k)*/
        if (k-j>MINGAP && compat(Enc->seq[j-1], Enc->seq[k-1])) {
          Enc->bp_left=-i;
          Enc->bp_right=-j;
          Enc->bp_left2=j;
          Enc->bp_right2=k;
          cnt += update_deepest(Enc, str, minim);
          /* in case useFirst is on and structure is found, end*/
          if (Enc->first && cnt > 0) return cnt;
        }
      }

      if (Enc->verbose_lvl>1)
        vrna_message_info(stderr, "%2d bracket %2d position, inner switch", brack_num+1, i);
      /* inner switch*/
      for (k=i+1; k<j; k++) {
        /* jump to end of the sub-bracketing*/
        if (pt[k]!=0 && pt[k]>k/*'('*/) {
            k=pt[k];
            continue;
        }

        /* left switch (i,j) to (k,j)*/
        if (j-k>MINGAP && compat(Enc->seq[k-1], Enc->seq[j-1])) {
          Enc->bp_left=-i;
          Enc->bp_right=-j;
          Enc->bp_left2=k;
          Enc->bp_right2=j;
          cnt += update_deepest(Enc, str, minim);
          /* in case useFirst is on and structure is found, end*/
          if (Enc->first && cnt > 0) return cnt;
        }

        /* right switch (i,j) to (i,k)*/
        if (k-i>MINGAP && compat(Enc->seq[i-1], Enc->seq[k-1])) {
          Enc->bp_left=-i;
          Enc->bp_right=-j;
          Enc->bp_left2=i;
          Enc->bp_right2=k;
          cnt += update_deepest(Enc, str, minim);
          /* in case useFirst is on and structure is found, end*/
          if (Enc->first && cnt > 0) return cnt;
        }
      } /* end inner switch for*/
      brack_num++;
    } /* end if (pt[i]=='(')*/
  } /* end for in switches*/
  return cnt;
}

/* move to deepest (or first) neighbour*/
PRIVATE int
move_set(Encoded *Enc, struct_en *str){

  /* count how many times called*/
  cnt_move++;

  /* count better neighbours*/
  int cnt = 0;

  /* deepest descent*/
  struct_en min;
  min.structure = allocopy(str->structure);
  min.energy = str->energy;
  Enc->current_en = str->energy;

  if (Enc->verbose_lvl>0) { fprintf(stderr, "  start of MS:\n  "); print_str(stderr, str->structure); fprintf(stderr, " %d\n\n", str->energy); }

  /* if using first dont do all of them*/
  bool end = false;
  /* insertions*/
  if (!end) cnt += insertions(Enc, str, &min);
  if (Enc->first && cnt>0) end = true;
  if (Enc->verbose_lvl>1) fprintf(stderr, "\n");

  /* deletions*/
  if (!end) cnt += deletions(Enc, str, &min);
  if (Enc->first && cnt>0) end = true;

  /* shifts (only if enabled + noLP disabled)*/
  if (!end && Enc->shift && !Enc->noLP) {
    cnt += shifts(Enc, str, &min);
    if (Enc->first && cnt>0) end = true;
  }

  /* if degeneracy occurs, solve it!*/
  if (!end && (Enc->end_unpr - Enc->begin_unpr)>0) {
    Enc->processed[Enc->end_pr] = str->structure;
    Enc->end_pr++;
    str->structure = Enc->unprocessed[Enc->begin_unpr];
    Enc->unprocessed[Enc->begin_unpr]=NULL;
    Enc->begin_unpr++;
    cnt += move_set(Enc, str);
  } else {
    /* write output to str*/
    copy_arr(str->structure, min.structure);
    str->energy = min.energy;
  }
  /* release minimal*/
  free(min.structure);

  /* resolve degeneracy in local minima*/
  if ((Enc->end_pr - Enc->begin_pr)>0) {
    Enc->processed[Enc->end_pr]=str->structure;
    Enc->end_pr++;

    int min = find_min(Enc->processed, Enc->begin_pr, Enc->end_pr);
    short *tmp = Enc->processed[min];
    Enc->processed[min] = Enc->processed[Enc->begin_pr];
    Enc->processed[Enc->begin_pr] = tmp;
    str->structure = Enc->processed[Enc->begin_pr];
    Enc->begin_pr++;
    free_degen(Enc);
  }

  if (Enc->verbose_lvl>1 && !(Enc->first)) { fprintf(stderr, "\n  end of MS:\n  "); print_str(stderr, str->structure); fprintf(stderr, " %d\n\n", str->energy); }

  return cnt;
}

PRIVATE void
construct_moves(Encoded *Enc, short *structure){

  /* generate all possible moves (less than n^2)*/
  Enc->num_moves = 0;
  int i;
  for (i=1; i<=structure[0]; i++) {
    if (structure[i]!=0) {
      if (structure[i]<i) continue;
      Enc->moves_from[Enc->num_moves]=-i;
      Enc->moves_to[Enc->num_moves]=-structure[i];
      Enc->num_moves++;
      /* fprintf(stderr, "add  d(%d, %d)\n", i, str.structure[i]); */
    } else {
      int j;
      for (j=i+1; j<=structure[0]; j++) {
        /* fprintf(stderr, "check (%d, %d)\n", i, j); */
        if (structure[j]==0) {
          if (try_insert_seq(Enc->seq,i,j)) {
            Enc->moves_from[Enc->num_moves]=i;
            Enc->moves_to[Enc->num_moves]=j;
            Enc->num_moves++;
            /* fprintf(stderr, "add  i(%d, %d)\n", i, j); */
            continue;
          }
        } else if (structure[j]>j) { /*  '(' */
          j = structure[j];
        } else break;
      }
    }
  }

  /* permute them */
  for (i=0; i<Enc->num_moves-1; i++) {
    int rnd = rand();
    rnd = rnd % (Enc->num_moves-i) + i;
    int swp;
    swp = Enc->moves_from[i];
    Enc->moves_from[i]=Enc->moves_from[rnd];
    Enc->moves_from[rnd]=swp;
    swp = Enc->moves_to[i];
    Enc->moves_to[i]=Enc->moves_to[rnd];
    Enc->moves_to[rnd]=swp;
  }
}

PRIVATE int
move_rset(Encoded *Enc, struct_en *str){

  /* count how many times called*/
  cnt_move++;

  /* count better neighbours*/
  int cnt = 0;

  /* deepest descent*/
  struct_en min;
  min.structure = allocopy(str->structure);
  min.energy = str->energy;
  Enc->current_en = str->energy;

  if (Enc->verbose_lvl>0) { fprintf(stderr, "  start of MR:\n  "); print_str(stderr, str->structure); fprintf(stderr, " %d\n\n", str->energy); }

  /*  construct and permute possible moves */
  construct_moves(Enc, str->structure);

  /* find first lower one*/
  int i;
  for (i=0; i<Enc->num_moves; i++) {
    Enc->bp_left = Enc->moves_from[i];
    Enc->bp_right = Enc->moves_to[i];
    cnt = update_deepest(Enc, str, &min);
    if (cnt) break;
  }

  /* if degeneracy occurs, solve it!*/
  if (!cnt && (Enc->end_unpr - Enc->begin_unpr)>0) {
    Enc->processed[Enc->end_pr] = str->structure;
    Enc->end_pr++;
    str->structure = Enc->unprocessed[Enc->begin_unpr];
    Enc->unprocessed[Enc->begin_unpr]=NULL;
    Enc->begin_unpr++;
    cnt += move_rset(Enc, str);
  } else {
    /* write output to str*/
    copy_arr(str->structure, min.structure);
    str->energy = min.energy;
  }
  /* release minimal*/
  free(min.structure);

  /* resolve degeneracy in local minima*/
  if ((Enc->end_pr - Enc->begin_pr)>0) {
    Enc->processed[Enc->end_pr]=str->structure;
    Enc->end_pr++;

    int min = find_min(Enc->processed, Enc->begin_pr, Enc->end_pr);
    short *tmp = Enc->processed[min];
    Enc->processed[min] = Enc->processed[Enc->begin_pr];
    Enc->processed[Enc->begin_pr] = tmp;
    str->structure = Enc->processed[Enc->begin_pr];
    Enc->begin_pr++;
    free_degen(Enc);
  }

  return cnt;
}

/*check if base is lone*/
PRIVATE int
lone_base(short *pt, int i){

  if (i<=0 || i>pt[0]) return 0;
  /* is not a base pair*/
  if (pt[i]==0) return 0;

  /* base is lone:*/
  if (i-1>0) {
    /* is base pair and is the same bracket*/
    if (pt[i-1]!=0 && ((pt[i-1]<pt[pt[i-1]]) == (pt[i]<pt[pt[i]]))) return 0;
  }

  if (i+1<=pt[0]) {
    if (pt[i+1]!=0 && ((pt[i-1]<pt[pt[i-1]]) == (pt[i]<pt[pt[i]]))) return 0;
  }

  return 1;
}

/* if the structure has lone pairs*/
PRIVATE int
find_lone_pair(short* str){

  int i;
  for(i=1; i<str[0]; i++) {
    if (str[i]==0) continue; /* '.'*/

    if (str[i]>str[str[i]]) {  /* '('*/
      if (i+1==str[0] || str[i+1]==0 || str[i+1]<str[str[i+1]]) {
        return i;
      } else while (i+1!=str[0] && str[i+1]!=0 && str[i+1]>str[str[i+1]]) i++;
    }

    if (str[i]<str[str[i]]) {  /* ')'*/
      if (i+1==str[0] || str[i+1]==0 || str[i+1]>str[str[i+1]]) {
        return i;
      } else while (i+1!=str[0] && str[i+1]!=0 && str[i+1]<str[str[i+1]]) i++;
    }
  }

  return -1;
}

PUBLIC int
move_standard(char *seq,
              char *struc,
              enum MOVE_TYPE type,
              int verbosity_level,
              int shifts,
              int noLP){

  make_pair_matrix();

  short int *s0 = encode_sequence(seq, 0);
  short int *s1 = encode_sequence(seq, 1);
  short int *str = vrna_ptable(struc);

  int energy = 0;
  switch (type){
  case GRADIENT: energy = move_gradient(seq, str, s0, s1, verbosity_level, shifts, noLP); break;
  case FIRST: energy = move_first(seq, str, s0, s1, verbosity_level, shifts, noLP); break;
  case ADAPTIVE: energy = move_adaptive(seq, str, s0, s1, verbosity_level); break;
  }

  int i=1;
  for (; i<=str[0]; i++) {
    if (str[i]==0) struc[i-1]='.';
    else if (str[i]>str[str[i]]) struc[i-1]='(';
      else struc[i-1]=')';
  }

  free(s0);
  free(s1);
  free(str);

  return energy;
}

PUBLIC int
move_gradient(char *string,
              short *ptable,
              short *s,
              short *s1,
              int verbosity_level,
              int shifts,
              int noLP){

  cnt_move = 0;

  Encoded enc;
  enc.seq = string;
  enc.s0 = s;
  enc.s1 = s1;

  /* moves*/
  enc.bp_left=0;
  enc.bp_right=0;
  enc.bp_left2=0;
  enc.bp_right2=0;

  /* options*/
  enc.noLP=noLP;
  enc.verbose_lvl=verbosity_level;
  enc.first=0;
  enc.shift=shifts;

  /* degeneracy*/
  enc.begin_unpr=0;
  enc.begin_pr=0;
  enc.end_unpr=0;
  enc.end_pr=0;
  enc.current_en=0;

  /*  function */
  enc.funct=NULL;

  int i;
  for (i=0; i<MAX_DEGEN; i++) enc.processed[i]=enc.unprocessed[i]=NULL;

  struct_en str;
  str.structure = allocopy(ptable);
  str.energy = energy_of_structure_pt(enc.seq, str.structure, enc.s0, enc.s1, 0);

  while (move_set(&enc, &str)!=0) {
    free_degen(&enc);
  }
  free_degen(&enc);

  copy_arr(ptable, str.structure);
  free(str.structure);

  return str.energy;
}

PUBLIC int
move_first( char *string,
            short *ptable,
            short *s,
            short *s1,
            int verbosity_level,
            int shifts,
            int noLP){

  cnt_move = 0;

  Encoded enc;
  enc.seq = string;
  enc.s0 = s;
  enc.s1 = s1;

  /* moves*/
  enc.bp_left=0;
  enc.bp_right=0;
  enc.bp_left2=0;
  enc.bp_right2=0;

  /* options*/
  enc.noLP=noLP;
  enc.verbose_lvl=verbosity_level;
  enc.first=1;
  enc.shift=shifts;

  /* degeneracy*/
  enc.begin_unpr=0;
  enc.begin_pr=0;
  enc.end_unpr=0;
  enc.end_pr=0;
  enc.current_en=0;

  /*  function */
  enc.funct=NULL;

  int i;
  for (i=0; i<MAX_DEGEN; i++) enc.processed[i]=enc.unprocessed[i]=NULL;

  struct_en str;
  str.structure = allocopy(ptable);
  str.energy = energy_of_structure_pt(enc.seq, str.structure, enc.s0, enc.s1, 0);

  while (move_set(&enc, &str)!=0) {
    free_degen(&enc);
  }
  free_degen(&enc);

  copy_arr(ptable, str.structure);
  free(str.structure);

  return str.energy;
}

PUBLIC int
move_adaptive(char *string,
              short *ptable,
              short *s,
              short *s1,
              int verbosity_level){

  srand(time(NULL));

  cnt_move = 0;

  Encoded enc;
  enc.seq = string;
  enc.s0 = s;
  enc.s1 = s1;

  /* moves*/
  enc.bp_left=0;
  enc.bp_right=0;
  enc.bp_left2=0;
  enc.bp_right2=0;

  /* options*/
  enc.noLP=0;
  enc.verbose_lvl=verbosity_level;
  enc.first=1;
  enc.shift=0;

  /* degeneracy*/
  enc.begin_unpr=0;
  enc.begin_pr=0;
  enc.end_unpr=0;
  enc.end_pr=0;
  enc.current_en=0;

  /*  function */
  enc.funct=NULL;

  /*  allocate memory for moves */
  enc.moves_from = (int*) vrna_alloc(ptable[0]*ptable[0]*sizeof(int));
  enc.moves_to = (int*) vrna_alloc(ptable[0]*ptable[0]*sizeof(int));

  int i;
  for (i=0; i<MAX_DEGEN; i++) enc.processed[i]=enc.unprocessed[i]=NULL;

  struct_en str;
  str.structure = allocopy(ptable);
  str.energy = energy_of_structure_pt(enc.seq, str.structure, enc.s0, enc.s1, 0);

  while (move_rset(&enc, &str)!=0) {
    free_degen(&enc);
  }
  free_degen(&enc);

  copy_arr(ptable, str.structure);
  free(str.structure);
  free(enc.moves_from);
  free(enc.moves_to);

  return str.energy;
}

PUBLIC int
browse_neighs(char *seq,
              char *struc,
              int verbosity_level,
              int shifts,
              int noLP,
              int (*funct) (struct_en*, struct_en*)){

  make_pair_matrix();

  short int *s0 = encode_sequence(seq, 0);
  short int *s1 = encode_sequence(seq, 1);
  short int *str = vrna_ptable(struc);

  int res = browse_neighs_pt(seq, str, s0, s1, verbosity_level, shifts, noLP, funct);

  free(s0);
  free(s1);
  free(str);

  return res;
}

PUBLIC int
browse_neighs_pt( char *string,
                  short *ptable,
                  short *s,
                  short *s1,
                  int verbosity_level,
                  int shifts,
                  int noLP,
                  int (*funct) (struct_en*, struct_en*)){

  cnt_move = 0;

  Encoded enc;
  enc.seq = string;
  enc.s0 = s;
  enc.s1 = s1;

  /* moves*/
  enc.bp_left=0;
  enc.bp_right=0;
  enc.bp_left2=0;
  enc.bp_right2=0;

  /* options*/
  enc.noLP=noLP;
  enc.verbose_lvl=verbosity_level;
  enc.first=1;
  enc.shift=shifts;

  /* degeneracy*/
  enc.begin_unpr=0;
  enc.begin_pr=0;
  enc.end_unpr=0;
  enc.end_pr=0;
  enc.current_en=0;

  /*  function */
  enc.funct=funct;

  int i;
  for (i=0; i<MAX_DEGEN; i++) enc.processed[i]=enc.unprocessed[i]=NULL;

  struct_en str;
  str.structure = allocopy(ptable);
  str.energy = energy_of_structure_pt(enc.seq, str.structure, enc.s0, enc.s1, 0);

  move_set(&enc, &str);
  free_degen(&enc);

  copy_arr(ptable, str.structure);
  free(str.structure);

  return str.energy;
}

/* printf*/
PUBLIC void
print_stren(FILE *out, struct_en *str) {
  print_str(out, str->structure);
  fprintf(out, " %6.2f\n", str->energy/100.0);
}

PUBLIC void
print_str(FILE *out, short *str) {
  int i;
  for (i=1; i<=str[0]; i++) {
    if (str[i]==0) fprintf(out, ".");
    else if (str[i]<i) fprintf(out, ")");
    else fprintf(out, "(");
  }
}


#ifdef TEST_MOVESET
/*  sample usage: */
int main() {
  char seq[20] = "ACCCCCCTCTGTAGGGGGA";
  char str[20] = ".((.(.........).)).";

  /*  move to the local minimum and display it */
  int energy = move_standard(seq, str, GRADIENT, 0, 0, 0);
  fprintf(stdout, "%s %6.2f\n\n", str, energy/100.0);

  /* now create an array of every structure in neighbourhood of str structure */
  struct_en *list = NULL;
  int list_length = 0;

  int get_list(struct_en *new_one, struct_en *old_one)
  {
    /*  enlarge the list */
    list_length++;
    list = (struct_en*) realloc(list, list_length*sizeof(struct_en));

    /*  copy the structure */
    list[list_length-1].energy = new_one->energy;
    list[list_length-1].structure = allocopy(new_one->structure);

    /*  we want to continue -> return 0 */
    return 0;
  }
  browse_neighs(seq, str, 0, 0, 0, get_list);

  /*  print them and free the memory: */
  int i;
  for (i=0; i<list_length; i++) {
    print_stren(stdout, &list[i]);
    free(list[i].structure);
  }
  free(list);

  return 0;
}

#endif