/* Last changed Time-stamp: <2007-10-30 14:06:22 htafer> */
/*                
           compute the duplex structure of two RNA strands,
                allowing only inter-strand base pairs.
         see cofold() for computing hybrid structures without
                             restriction.
                             Ivo Hofacker
                          Vienna RNA package

*/


/*
  library containing the function used in rnaplex
  the program rnaplex uses the following function
  Lduplexfold: finds high scoring segments
  it stores the end-position of these segments in an array 
  and call then for each of these positions the duplexfold function 
  which allows to make backtracking for each of the high scoring position
  It allows to find suboptimal partially overlapping (depends on a a parameter) 
  duplexes between a long RNA and a shorter one.
  Contrarly to RNAduplex, the energy model is not in E~log(N), 
  where N is the length of an interial loop but used an affine model, 
  where the extension and begin parameter are fitted to the energy
  parameter used by RNAduplex. This allows to check for duplex between a short RNA(20nt) 
  and a long one at the speed of 1Mnt/s. At this speed the whole genome (3Gnt) can be analyzed for one siRNA 
  in about 50 minutes.
  The algorithm is based on an idea by Durbin and Eddy:when the alginment reach a value larger than a
  given threshold this value is stored in an array. When the alignment score goes 
  then under this threshold, the alignemnent begin from this value, in that way the backtracking allow us 
  to find all non-overlapping high-scoring segments. 
  For more information check "durbin, biological sequence analysis"
*/

#include <config.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <ctype.h>
#include <string.h>
#include "utils.h"
#include "energy_par.h"
#include "fold_vars.h"
#include "fold.h"
#include "pair_mat.h"
#include "params.h"
#include "loop_energies.h"
#include "plex.h"
#include "ali_plex.h"

/*@unused@*/
static char rcsid[] UNUSED = "$Id: plex.c,v 1.14 2007/06/12 12:50:16 htafer Exp $";


#define PUBLIC
#define PRIVATE static

#define STACK_BULGE1  1   /* stacking energies for bulges of size 1 */
#define NEW_NINIO     1   /* new asymetry penalty */
#define ARRAY 32          /*array size*/
#define UNIT 100
#define MINPSCORE -2 * UNIT
/**
*** Due to the great similarity between functions,
*** more annotation can be found in plex.c
**/

PRIVATE short *encode_seq(const char *seq);
PRIVATE void  update_dfold_params(void);
/**
*** aliduplexfold(_XS)/alibacktrack(_XS) computes duplex interaction with standard energy and considers extension_cost 
*** alifind_max(_XS)/aliplot_max(_XS) find suboptimals and MFE
**/
PRIVATE duplexT aliduplexfold(const char *s1[], const char *s2[], const int extension_cost);
PRIVATE char *    alibacktrack(int i, int j, const short *s1[], const short *s2[], const int extension_cost);
PRIVATE void      alifind_max(const int *position, const int *position_j,const int delta, const int threshold, 
                           const int alignment_length, const char *s1[], const char *s2[], const int extension_cost, const int fast);
PRIVATE void      aliplot_max(const int max, const int max_pos, const int max_pos_j, 
                           const int alignment_length, const char *s1[], const char *s2[], const int extension_cost, const int fast);
PRIVATE duplexT aliduplexfold_XS(const char *s1[], const char *s2[],const int **access_s1, 
                                const int **access_s2, const int i_pos, const int j_pos, const int threshold);
PRIVATE char *    alibacktrack_XS(int i, int j, const short *s1[], const short *s2[], const int** access_s1, const int** access_s2);
PRIVATE void  alifind_max_XS(const int *position, const int *position_j, 
                                 const int delta,  const int threshold, const int alignment_length,
                                 const char* s1[], const char* s2[], 
                                 const int **access_s1, const int **access_s2, const int fast);
PRIVATE void aliplot_max_XS(const int max, const int max_pos, const int max_pos_j,
                            const int alignment_length, const char *s1[], const char* s2[], 
                            const int **access_s1, const int **access_s2, const int fast);

/**
*** computes covariance score
**/

PRIVATE int covscore(const int *types, int n_seq);

extern double cv_fact; /* from alifold.c, default 1 */
extern double nc_fact;


/*@unused@*/

#define MAXSECTORS      500     /* dimension for a backtrack array */
#define LOCALITY        0.      /* locality parameter for base-pairs */

PRIVATE paramT *P = NULL;
PRIVATE int   **c = NULL;
PRIVATE int  **lc = NULL, **lin = NULL, **lbx = NULL, **lby = NULL,**linx = NULL, **liny = NULL;   
                                             



PRIVATE int   n1,n2;
PRIVATE int n3, n4; 
PRIVATE int delay_free=0;


/*-----------------------------------------------------------------------duplexfold_XS---------------------------------------------------------------------------*/



/*----------------------------------------------ALIDUPLEXFOLD-----------------------------------------------------------------------------------------------------------*/
PRIVATE duplexT aliduplexfold(const char *s1[], const char *s2[], const int extension_cost) {
  int i, j, s, n_seq, l1, Emin=INF, i_min=0, j_min=0;
  char *struc;
  duplexT mfe;
  short **S1, **S2;
  int *type;
  n3 = (int) strlen(s1[0]);
  n4 = (int) strlen(s2[0]);
  for (s=0; s1[s]!=NULL; s++);
  n_seq = s;
  for (s=0; s2[s]!=NULL; s++);
  if (n_seq != s) nrerror("unequal number of sequences in aliduplexfold()\n");

  if ((!P) || (fabs(P->temperature - temperature)>1e-6)) {
    update_fold_params();  if(P) free(P); P = scale_parameters();
    make_pair_matrix();
  }
  
  c = (int **) space(sizeof(int *) * (n3+1));
  for (i=1; i<=n3; i++) c[i] = (int *) space(sizeof(int) * (n4+1));
  
  S1 = (short **) space((n_seq+1)*sizeof(short *));
  S2 = (short **) space((n_seq+1)*sizeof(short *));
  for (s=0; s<n_seq; s++) {
    if (strlen(s1[s]) != n3) nrerror("uneqal seqence lengths");
    if (strlen(s2[s]) != n4) nrerror("uneqal seqence lengths");
    S1[s] = encode_seq(s1[s]);
    S2[s] = encode_seq(s2[s]);
  }
  type = (int *) space(n_seq*sizeof(int));

  for (i=1; i<=n3; i++) {
    for (j=n4; j>0; j--) {
      int k,l,E,psc;
      for (s=0; s<n_seq; s++) {
        type[s] = pair[S1[s][i]][S2[s][j]];
      }
      psc = covscore(type, n_seq);
      for (s=0; s<n_seq; s++) if (type[s]==0) type[s]=7;
      c[i][j] = (psc>=MINPSCORE) ? (n_seq*(P->DuplexInit + 2*extension_cost)) : INF;
      if (psc<MINPSCORE) continue;
      for (s=0; s<n_seq; s++) {
        c[i][j] += E_ExtLoop(type[s], (i>1) ? S1[s][i-1] : -1, (j<n4) ? S2[s][j+1] : -1, P) + 2*extension_cost;
      }
      for (k=i-1; k>0 && k>i-MAXLOOP-2; k--) {
        for (l=j+1; l<=n4; l++) {
          int type2;
          if (i-k+l-j-2>MAXLOOP) break;
          if (c[k][l]>INF/2) continue;
          for (E=s=0; s<n_seq; s++) { 
            type2 = pair[S1[s][k]][S2[s][l]];
            if (type2==0) type2=7;
            E += E_IntLoop(i-k-1, l-j-1, type2, rtype[type[s]],
                           S1[s][k+1], S2[s][l-1], S1[s][i-1], S2[s][j+1],P) + (i-k+l-j)*extension_cost;
          }
          c[i][j] = MIN2(c[i][j], c[k][l]+E);
        }
      }
      c[i][j] -= psc;
      E = c[i][j]; 
      for (s=0; s<n_seq; s++) {
        E += E_ExtLoop(rtype[type[s]], (j>1) ? S2[s][j-1] : -1, (i<n3) ? S1[s][i+1] : -1, P) +2*extension_cost;
      }
      if (E<Emin) {
        Emin=E; i_min=i; j_min=j;
      } 
    }
  }
  struc = alibacktrack(i_min, j_min, (const short int**) S1, (const short int**) S2 , extension_cost);
  if (i_min<n3) i_min++;
  if (j_min>1 ) j_min--;
  l1 = strchr(struc, '&')-struc;
  int size; 
  size=strlen(struc)-1;
  Emin-=size * n_seq * extension_cost;
  mfe.i = i_min;
  mfe.j = j_min;
  mfe.energy = (float) (Emin/(100.*n_seq));
  mfe.structure = struc;
  if (!delay_free) {
    for (i=1; i<=n3; i++) free(c[i]);
    free(c);
  }
  for (s=0; s<n_seq; s++) {
    free(S1[s]); free(S2[s]);
  }
  free(S1); free(S2); free(type);
  return mfe;
}


PRIVATE char *alibacktrack(int i, int j, const short *S1[], const short *S2[], const int extension_cost) {
  /* backtrack structure going backwards from i, and forwards from j 
     return structure in bracket notation with & as separator */
  int k, l, *type, type2, E, traced, i0, j0, s, n_seq;
  char *st1, *st2, *struc;
  
  n3 = (int) S1[0][0];
  n4 = (int) S2[0][0];

  for (s=0; S1[s]!=NULL; s++);
  n_seq = s;
  for (s=0; S2[s]!=NULL; s++);
  if (n_seq != s) nrerror("unequal number of sequences in alibacktrack()\n");

  st1 = (char *) space(sizeof(char)*(n3+1));
  st2 = (char *) space(sizeof(char)*(n4+1));
  type = (int *) space(n_seq*sizeof(int));

  i0=MIN2(i+1,n3); j0=MAX2(j-1,1);

  while (i>0 && j<=n4) {
    int psc;
    E = c[i][j]; traced=0;
    st1[i-1] = '(';
    st2[j-1] = ')'; 
    for (s=0; s<n_seq; s++) {
      type[s] = pair[S1[s][i]][S2[s][j]];
    }
    psc = covscore(type, n_seq);
    for (s=0; s<n_seq; s++) if (type[s]==0) type[s] = 7;
    E += psc;
    for (k=i-1; k>0 && k>i-MAXLOOP-2; k--) {
      for (l=j+1; l<=n4; l++) {
        int LE;
        if (i-k+l-j-2>MAXLOOP) break;
        if (c[k][l]>INF/2) continue;
        for (s=LE=0; s<n_seq; s++) {
          type2 = pair[S1[s][k]][S2[s][l]];
          if (type2==0) type2=7;
          LE += E_IntLoop(i-k-1, l-j-1, type2, rtype[type[s]],
                          S1[s][k+1], S2[s][l-1], S1[s][i-1], S2[s][j+1],P)+(i-k+l-j)*extension_cost;
        }
        if (E == c[k][l]+LE) {
          traced=1; 
          i=k; j=l;
          break;
        }
      }
      if (traced) break;
    }
    if (!traced) { 
      for (s=0; s<n_seq; s++) {
        E -= E_ExtLoop(type[s], (i>1) ? S1[s][i-1] : -1, (j<n4) ? S2[s][j+1] : -1, P) + 2*extension_cost;
      }
      if (E != n_seq*P->DuplexInit + n_seq*2*extension_cost) {
        nrerror("backtrack failed in aliduplex");
      } else break;
    }
  }
  if (i>1)  i--;
  if (j<n4) j++;
  
  struc = (char *) space(i0-i+1+j-j0+1+2);
  for (k=MAX2(i,1); k<=i0; k++) if (!st1[k-1]) st1[k-1] = '.';
  for (k=j0; k<=j; k++) if (!st2[k-1]) st2[k-1] = '.';
  strcpy(struc, st1+MAX2(i-1,0)); strcat(struc, "&"); 
  strcat(struc, st2+j0-1);
  
  /* printf("%s %3d,%-3d : %3d,%-3d\n", struc, i,i0,j0,j);  */
  free(st1); free(st2); free(type);

  return struc;
}

duplexT** aliLduplexfold(const char *s1[], const char *s2[], const int threshold, const int extension_cost, const int alignment_length, const int delta, const int fast,const int il_a, const int il_b, const int b_a, const int b_b)
{
  short **S1, **S2;
  int *type, type2;
  int i, j,s,n_seq;
  s=0;
  int bopen=b_b;
  int bext=b_a+extension_cost;
  int iopen=il_b;
  int iext_s=2*(il_a+extension_cost);//iext_s 2 nt nucleotide extension of interior loop, on i and j side
  int iext_ass=50+il_a+extension_cost;//iext_ass assymetric extension of interior loop, either on i or on j side.
  int min_colonne=INF; //enthaelt das maximum einer kolonne
  int i_length;
  int max_pos;//get position of the best hit
  int max_pos_j;
  int temp;
  int min_j_colonne;
  int max=INF;
  //FOLLOWING NEXT 4 LINE DEFINES AN ARRAY CONTAINING POSITION OF THE SUBOPT IN S1
  int *position; //contains the position of the hits with energy > E
  int *position_j;
  
  
  n1 = (int) strlen(s1[0]);
  n2 = (int) strlen(s2[0]);
  for (s=0; s1[s]; s++);
  n_seq = s;
  for (s=0; s2[s]; s++);
  if (n_seq != s) nrerror("unequal number of sequences in aliduplexfold()\n");
  
  position = (int *) space ((delta+(n1)+4+delta) * sizeof(int));
  position_j= (int *) space ((delta+(n1)+4+delta) * sizeof(int));
  
  if ((!P) || (fabs(P->temperature - temperature)>1e-6)){
    update_dfold_params();
  }
  
  lc   = (int**) space (sizeof(int *) * 5);
  lin  = (int**) space (sizeof(int *) * 5);
  lbx  = (int**) space (sizeof(int *) * 5);
  lby  = (int**) space (sizeof(int *) * 5);
  linx = (int**) space (sizeof(int *) * 5);
  liny = (int**) space (sizeof(int *) * 5);

  for (i=0; i<=4; i++){
    lc[i]  = (int *) space(sizeof(int) * (n2+5));  
    lin[i] = (int *) space(sizeof(int) * (n2+5));  
    lbx[i] = (int *) space(sizeof(int) * (n2+5));  
    lby[i] = (int *) space(sizeof(int) * (n2+5));  
    linx[i]= (int *) space(sizeof(int) * (n2+5));  
    liny[i]= (int *) space(sizeof(int) * (n2+5));  
  }

  
  S1 = (short **) space((n_seq+1)*sizeof(short *));
  S2 = (short **) space((n_seq+1)*sizeof(short *));
  for (s=0; s<n_seq; s++) {
    if (strlen(s1[s]) != n1) nrerror("uneqal seqence lengths");
    if (strlen(s2[s]) != n2) nrerror("uneqal seqence lengths");
    S1[s] = encode_seq(s1[s]);
    S2[s] = encode_seq(s2[s]);    
  }
  type = (int *) space(n_seq*sizeof(int));
  /**
  *** array initialization
  **/
  for(j=n2;j>=0;j--) {
    lbx[0][j]=lbx[1][j]=lbx[2][j]=lbx[3][j]    = lbx[4][j] =INF;
    lin[0][j]=lin[1][j]=lin[2][j]=lin[3][j]    = lin[4][j] =INF;
    lc[0][j] =lc[1][j] =lc[2][j] = lc[3][j]    =  lc[4][j] =INF;
    lby[0][j]=lby[1][j]=lby[2][j]=lby[3][j]    = lby[4][j] =INF;
    liny[0][j]=liny[1][j]=liny[2][j]=liny[3][j]=liny[4][j]=INF;
    linx[0][j]=linx[1][j]=linx[2][j]=linx[3][j]=linx[4][j]=INF;
  }
  i=10;
  i_length= n1 - 9 ;
  while(i < i_length) {
    int idx=i%5;
    int idx_1=(i-1)%5;
    int idx_2=(i-2)%5;
    int idx_3=(i-3)%5;
    int idx_4=(i-4)%5;
    j=n2-9;
    while (9 < --j) {
      int psc;
      for (s=0; s<n_seq; s++) {
        type[s] = pair[S1[s][i]][S2[s][j]];
      }
      psc = covscore(type, n_seq);
      for (s=0; s<n_seq; s++) if (type[s]==0) type[s]=7;
      lc[idx][j] = (psc>=MINPSCORE) ? (n_seq*P->DuplexInit + 2*n_seq*extension_cost) : INF;
      /**
      *** Update matrix. It is the average over all sequence of a given structure element
      *** c_stack -> stacking of c
      *** c_10, c01 -> stack from bulge
      *** c_nm -> arrives in stack from nxm loop
      *** c_in -> arrives in stack from interior loop
      *** c_bx -> arrives in stack from large bulge on target
      *** c_by -> arrives in stack from large bulge on query
      *** 
      **/
      int c_stack, c_10, c_01, c_11, c_22, c_21, c_12, c_23, c_32, c_in, c_in2x, c_in2y, c_bx, c_by, c_inx, c_iny;  //matrix c
      int in, in_x, in_y, in_xy; // in begin, in_x assymetric, in_y assymetric, in_xy symetric;
      int inx, inx_x;
      int iny, iny_y;
      int bx, bx_x; 
      int by, by_y; 
      in=lc[idx_1][j+1]; in_x=lin[idx_1][j]; in_y=lin[idx][j+1]; in_xy=lin[idx_1][j+1];
      inx=lc[idx_1][j+1]; inx_x=linx[idx_1][j];
      iny=lc[idx_1][j+1]; iny_y=liny[idx][j+1];
      bx=lc[idx_1][j]; bx_x=lbx[idx_1][j]; 
      by=lc[idx][j+1]; by_y=lby[idx][j+1];
      c_stack=lc[idx_1][j+1]; c_01=lc[idx_1][j+2];c_10=lc[idx_2][j+1]; 
      c_12=lc[idx_2][j+3];c_21=lc[idx_3][j+2];c_11=lc[idx_2][j+2];
      c_22=lc[idx_3][j+3];c_32=lc[idx_4][j+3];c_23=lc[idx_3][j+4];
      c_in=lin[idx_3][j+3];c_in2x=lin[idx_4][j+2];c_in2y=lin[idx_2][j+4];
      c_inx=linx[idx_3][j+1]; c_iny=liny[idx_1][j+3];
      c_bx=lbx[idx_2][j+1];c_by=lby[idx_1][j+2];
      for (s=0; s<n_seq; s++) {
        type2 = pair[S2[s][j+1]][S1[s][i-1]];
        in   +=P->mismatchI[type2][S2[s][j]][S1[s][i]]+iopen+iext_s;
        in_x +=iext_ass;
        in_y +=iext_ass;
        in_xy+=iext_s;
        inx  +=P->mismatch1nI[type2][S2[s][j]][S1[s][i]]+iopen+iext_s;
        inx_x+=iext_ass;
        iny  +=P->mismatch1nI[type2][S2[s][j]][S1[s][i]]+iopen+iext_s;
        iny_y+=iext_ass;
        type2=pair[S2[s][j]][S1[s][i-1]];   
        bx   +=bopen+bext+(type2>2?P->TerminalAU:0);
        bx_x +=bext;
        type2=pair[S2[s][j+1]][S1[s][i]];
        by   +=bopen+bext+(type2>2?P->TerminalAU:0);
        by_y +=bext;
      }
      lin [idx][j]=MIN2(in, MIN2(in_x, MIN2(in_y, in_xy)));       
      linx[idx][j]=MIN2(inx_x, inx);
      liny[idx][j]=MIN2(iny_y, iny);
      lby[idx][j] =MIN2(by, by_y);
      lbx[idx][j] =MIN2(bx, bx_x);
      
      if (psc<MINPSCORE) continue;  
      for (s=0; s<n_seq; s++) {
        lc[idx][j]+=E_ExtLoop(type[s], S1[s][i-1],S2[s][j+1], P) + 2*extension_cost;
      }
      for (s=0; s<n_seq; s++) {
        type2=pair[S1[s][i-1]][S2[s][j+1]];if (type2==0) type2=7;
        c_stack+=E_IntLoop(0,0,type2, rtype[type[s]],S1[s][i], S2[s][j], S1[s][i-1], S2[s][j+1], P)+2*extension_cost;
        type2=pair[S1[s][i-1]][S2[s][j+2]];if (type2==0) type2=7;
        c_01   +=E_IntLoop(0,1,type2, rtype[type[s]],S1[s][i], S2[s][j+1], S1[s][i-1], S2[s][j+1], P)+3*extension_cost;
        type2=pair[S1[s][i-2]][S2[s][j+1]]; if (type2==0) type2=7;
        c_10   +=E_IntLoop(1,0,type2, rtype[type[s]],S1[s][i-1], S2[s][j], S1[s][i-1], S2[s][j+1], P)+3*extension_cost;
        type2=pair[S1[s][i-2]][S2[s][j+2]]; if (type2==0) type2=7;
        c_11   +=E_IntLoop(1,1,type2, rtype[type[s]],S1[s][i-1], S2[s][j+1], S1[s][i-1], S2[s][j+1], P)+4*extension_cost;
        type2 = pair[S1[s][i-3]][S2[s][j+3]];if (type2==0) type2=7;
        c_22   +=E_IntLoop(2,2,type2, rtype[type[s]],S1[s][i-2], S2[s][j+2], S1[s][i-1], S2[s][j+1], P)+6*extension_cost;
        type2 = pair[S1[s][i-3]][S2[s][j+2]];if (type2==0) type2=7;
        c_21   +=E_IntLoop(2,1,type2, rtype[type[s]],S1[s][i-2], S2[s][j+1], S1[s][i-1], S2[s][j+1], P)+5*extension_cost;
        type2 = pair[S1[s][i-2]][S2[s][j+3]];if (type2==0) type2=7;
        c_12   +=E_IntLoop(1,2,type2, rtype[type[s]],S1[s][i-1], S2[s][j+2], S1[s][i-1], S2[s][j+1], P)+5*extension_cost;
        type2 = pair[S1[s][i-4]][S2[s][j+3]];if (type2==0) type2=7;
        c_32   +=E_IntLoop(3,2,type2, rtype[type[s]],S1[s][i-3], S2[s][j+2], S1[s][i-1], S2[s][j+1], P)+7*extension_cost;
        type2 = pair[S1[s][i-3]][S2[s][j+4]];if (type2==0) type2=7;
        c_23   +=E_IntLoop(2,3,type2, rtype[type[s]],S1[s][i-2], S2[s][j+3], S1[s][i-1], S2[s][j+1], P)+7*extension_cost;
        c_in   +=P->mismatchI[rtype[type[s]]][S1[s][i-1]][S2[s][j+1]]+2*extension_cost+2*iext_s;
        c_in2x +=P->mismatchI[rtype[type[s]]][S1[s][i-1]][S2[s][j+1]]+iext_s+2*iext_ass+2*extension_cost;
        c_in2y +=P->mismatchI[rtype[type[s]]][S1[s][i-1]][S2[s][j+1]]+iext_s+2*iext_ass+2*extension_cost;
        c_inx  +=P->mismatch1nI[rtype[type[s]]][S1[s][i-1]][S2[s][j+1]]+iext_ass+iext_ass+2*extension_cost;
        c_iny  +=P->mismatch1nI[rtype[type[s]]][S1[s][i-1]][S2[s][j+1]]+iext_ass+iext_ass+2*extension_cost;
        int bAU;
        bAU=(type[s]>2?P->TerminalAU:0);
        c_bx   +=2*extension_cost+bext+bAU;
        c_by   +=2*extension_cost+bext+bAU;
      }
      lc[idx][j] =MIN2(lc[idx][j], 
                       MIN2(c_stack, 
                            MIN2(c_10, 
                                 MIN2(c_01, 
                                      MIN2(c_11, 
                                           MIN2(c_21, 
                                                MIN2(c_12, 
                                                     MIN2(c_22,
                                                          MIN2(c_23,
                                                               MIN2(c_32,
                                                                    MIN2(c_bx, 
                                                                         MIN2(c_by, 
                                                                              MIN2(c_in,
                                                                                   MIN2(c_in2x,
                                                                                        MIN2(c_in2y,
                                                                                             MIN2(c_inx,c_iny)
                                                                                             )
                                                                                        )
                                                                                   )
                                                                              )
                                                                         )
                                                                    )
                                                               )
                                                          )
                                                     )
                                                )
                                           )
                                      )
                                 )
                            )
                       );
      lc[idx][j]-=psc;
      temp=lc[idx][j];
      for (s=0; s<n_seq; s++) {
        temp+=E_ExtLoop(rtype[type[s]], S2[s][j-1],S1[s][i+1],P)+2*extension_cost;
      }
      if(min_colonne > temp){
        min_colonne=temp;
        min_j_colonne=j;
      }
    }
    if(max>=min_colonne){
            max=min_colonne;
            max_pos=i;
        max_pos_j=min_j_colonne;
    }
    position[i+delta]=min_colonne;min_colonne=INF;
    position_j[i+delta]=min_j_colonne;
    i++;
  }
  //printf("MAX:%d ",max);
  for (s=0; s<n_seq; s++) {free(S1[s]);free(S2[s]);}
  free(S1); free(S2); 
  if(max<threshold){
    alifind_max(position, position_j, delta, threshold, alignment_length, s1, s2, extension_cost, fast);
  }
  aliplot_max(max, max_pos, max_pos_j,alignment_length, s1, s2, extension_cost,fast);
  for (i=0; i<=4; i++) {free(lc[i]);free(lin[i]);free(lbx[i]);free(lby[i]);free(linx[i]);free(liny[i]);}
  //free(lc[0]);free(lin[0]);free(lbx[0]);free(lby[0]);free(linx[0]);free(liny[0]);
  free(lc);free(lin);free(lbx);free(lby);free(linx);free(liny);
  free(position);
  free(position_j);
  free(type);
  return NULL;
}


PRIVATE void alifind_max(const int *position, const int *position_j,
                         const int delta, const int threshold, const int alignment_length, 
                         const char *s1[], const char *s2[], 
                         const int extension_cost, const int fast){
  int n_seq=0;
  for (n_seq=0; s1[n_seq]!=NULL; n_seq++);
  int pos=n1-9;
  if(fast==1){
    while(10<pos--){
      int temp_min=0;                                                               
      if(position[pos+delta]<(threshold)){                                          
        int search_range;                                                           
        search_range=delta+1;                                                       
        while(--search_range){                                                      
          if(position[pos+delta-search_range]<=position[pos+delta-temp_min]){                
            temp_min=search_range;                                                  
          }                                                                         
        }                                                                           
        pos-=temp_min;                                                              
        int max_pos_j;                                                              
        max_pos_j=position_j[pos+delta];
        int max;
        max=position[pos+delta];
        printf("target upper bound %d: query lower bound %d  (%5.2f) \n", pos-10, max_pos_j-10, ((double)max)/(n_seq*100));
        pos=MAX2(10,pos-delta);
      }
    }
  }
  else{
    pos=n1-9;
    while(pos-- > 10) {
      //printf("delta %d position:%d value:%d\n", delta, pos, position[pos]);
      int temp_min=0;
      if(position[pos+delta]<(threshold)){
        int search_range;
        search_range=delta+1;
        while(--search_range){
          if(position[pos+delta-search_range]<=position[pos+delta-temp_min]){
            temp_min=search_range;
          }
        }
        pos-=temp_min;
        int max_pos_j;
        max_pos_j=position_j[pos+delta];
        //printf("%d %d %d\n", pos, max_pos_j,position[pos+delta]);
        int begin_t=MAX2(11, pos-alignment_length+1);
        int end_t  =MIN2(n1-10, pos+1);
        int begin_q=MAX2(11, max_pos_j-1);
        int end_q  =MIN2(n2-10, max_pos_j+alignment_length-1);
        char **s3, **s4;
        s3 = (char**) space(sizeof(char*)*(n_seq+1));
        s4 = (char**) space(sizeof(char*)*(n_seq+1));
        int i;
        for(i=0; i<n_seq; i++){
          s3[i] = (char*) space(sizeof(char)*(end_t-begin_t+2));
          s4[i] = (char*) space(sizeof(char)*(end_q-begin_q+2));
          strncpy(s3[i], (s1[i]+begin_t-1), end_t - begin_t +1);
          strncpy(s4[i], (s2[i]+begin_q-1), end_q - begin_q +1);
          s3[i][end_t - begin_t +1]='\0';
          s4[i][end_q - begin_q +1]='\0';
        }
        duplexT test;
        test = aliduplexfold((const char**)s3, (const char**)s4, extension_cost);
        //printf("test %d threshold %d",test.energy*100,(threshold/n_seq));
        if(test.energy * 100 < (int) (threshold/n_seq)){
          int l1=strchr(test.structure, '&')-test.structure;
                   printf("%s %3d,%-3d : %3d,%-3d (%5.2f)\n", test.structure, 
                 begin_t -10 +test.i-l1, 
                 begin_t -10 +test.i-1, 
                 begin_q -10 + test.j - 1, 
                 begin_q-11 + test.j +strlen(test.structure) -l1 -2 , test.energy);
          pos=MAX2(10,pos-delta);
          
        }
        for(i=0;i<n_seq;i++){
          free(s3[i]);free(s4[i]);
        }
        free(s3);free(s4);
        free(test.structure);
      }
    }
  }
}
PRIVATE void aliplot_max(const int max, const int max_pos, const int max_pos_j, const int alignment_length, const char *s1[], const char *s2[], const int extension_cost, const int fast)
{
  int n_seq;
  for (n_seq=0; !(s1[n_seq]==NULL); n_seq++);
  n1 = strlen(s1[0]); //get length of alignment
  n2 = strlen(s2[0]); //get length of alignment
  if(fast==1){
    printf("target upper bound %d: query lower bound %d (%5.2f)\n", 
           max_pos-10, max_pos_j-10, (double) ((double)max)/(100*n_seq));
  }
  else{
    int begin_t=MAX2(11, max_pos-alignment_length+1);
    int end_t  =MIN2(n1-10, max_pos+1);
    int begin_q=MAX2(11, max_pos_j-1);
    int end_q  =MIN2(n2-10, max_pos_j+alignment_length-2);
    char **s3, **s4;
    s3 = (char**) space(sizeof(char*)*(n_seq+1));
    s4 = (char**) space(sizeof(char*)*(n_seq+1));
    int i;
    for(i=0; i<n_seq; i++){
      s3[i] = (char*) space(sizeof(char)*(end_t-begin_t+2));
      s4[i] = (char*) space(sizeof(char)*(end_q-begin_q+2));
      strncpy(s3[i], (s1[i]+begin_t-1), end_t - begin_t +1);
      strncpy(s4[i], (s2[i]+begin_q-1), end_q - begin_q +1);
      s3[i][end_t - begin_t +1]='\0';
      s4[i][end_q - begin_q +1]='\0';
    }
    duplexT test;
    s3[n_seq]=s4[n_seq]=NULL;
    test = aliduplexfold((const char**) s3,(const char**) s4, extension_cost);
    int l1=strchr(test.structure, '&')-test.structure;
    printf("%s %3d,%-3d : %3d,%-3d (%5.2f)\n", 
           test.structure, 
           begin_t -10 +test.i-l1, 
           begin_t -10 +test.i-1, 
           begin_q-10 + test.j - 1, 
           begin_q -11 + test.j +strlen(test.structure) - l1 - 2, 
           test.energy);
    for(i=0; i<n_seq ; i++){
      free(s3[i]);free(s4[i]);
    }
    free(s3);free(s4);
    free(test.structure);  
  }
}

PRIVATE duplexT aliduplexfold_XS(const char *s1[], const char *s2[],
                         const int **access_s1, const int **access_s2, 
                         const int i_pos, const int j_pos, const int threshold){
  int i,j,s,p,q, Emin=INF, l_min=0, k_min=0;
  char *struc;
  short **S1,**S2;
  int *type,*type2;
  struc=NULL;
  duplexT mfe;
  int n_seq;
  n3 = (int) strlen(s1[0]);
  n4 = (int) strlen(s2[0]);
  for (s=0; s1[s]!=NULL; s++);
  n_seq = s;
  for (s=0; s2[s]!=NULL; s++);
  //printf("%d \n",i_pos);
  if ((!P) || (fabs(P->temperature - temperature)>1e-6)) {
    update_fold_params();  if(P) free(P); P = scale_parameters();
    make_pair_matrix();
  }
  c = (int **) space(sizeof(int *) * (n3+1));
  for (i=0; i<=n3; i++) c[i] = (int *) space(sizeof(int) * (n4+1));
  for (i=0; i<=n3; i++){
    for(j=0;j<=n4;j++){
      c[i][j]=INF;
    }
  }
  S1 = (short **) space((n_seq+1)*sizeof(short *));
  S2 = (short **) space((n_seq+1)*sizeof(short *));
  for (s=0; s<n_seq; s++) {
    if (strlen(s1[s]) != n3) nrerror("uneqal seqence lengths");
    if (strlen(s2[s]) != n4) nrerror("uneqal seqence lengths");
    S1[s] = encode_seq(s1[s]);
    S2[s] = encode_seq(s2[s]);
  }
  type =  (int *) space(n_seq*sizeof(int));
  type2 = (int *) space(n_seq*sizeof(int));
  int type3, E, k,l;
  i=n3-1; j=2;
  for (s=0; s<n_seq; s++) {
    type[s] = pair[S1[s][i]][S2[s][j]];
  }
  c[i][j] = n_seq*P->DuplexInit - covscore(type,n_seq); 
  for (s=0; s<n_seq; s++) if (type[s]==0) type[s]=7; 
  for (s=0; s<n_seq; s++) {
    c[i][j]+=E_ExtLoop(rtype[type[s]], S2[s][j-1] ,S1[s][i+1],P);
  }
  for (k=i-1; k>0 ; k--) {
    c[k+1][0]=INF;
    for (l=j+1; l<=n4; l++) {
      c[k][l]=INF;
      int psc2;
      for(s=0;s<n_seq;s++){
        type2[s] = pair[S1[s][k]][S2[s][l]];
      }
      psc2=covscore(type2, n_seq);
      if (psc2<MINPSCORE) continue;
      for (s=0; s<n_seq; s++) if (type2[s]==0) type2[s]=7;
      for (p=k+1; p< n3 && p<k+MAXLOOP-1; p++) {                  
        for (q = l-1; q > 1; q--) {
          if (p-k+l-q-2>MAXLOOP) break;
          for(E=s=0;s<n_seq;s++){
            type3=pair[S1[s][p]][S2[s][q]];
            if(type3==0) type3=7;
            E += E_IntLoop(p-k-1, l-q-1, type2[s], rtype[type3],
                           S1[s][k+1], S2[s][l-1], S1[s][p-1], S2[s][q+1],P);
          }
          c[k][l] = MIN2(c[k][l], c[p][q]+E);
        }
      }
      c[k][l]-=psc2;
      E = c[k][l];
      E+=n_seq*(access_s1[i-k+1][i_pos]+access_s2[l-1][j_pos+(l-1)-1]);
      for (s=0; s<n_seq; s++) {
        E+=E_ExtLoop(type2[s], (k>1) ? S1[s][k-1] : -1, (l<n4) ? S2[s][l+1] : -1, P);
      }
      if (E<Emin) {
        Emin=E; k_min=k; l_min=l;
      } 
    }
  }
  if(Emin > threshold-1){ 
    mfe.structure=NULL; 
    mfe.energy=INF;
    for (i=0; i<=n3; i++) free(c[i]);
    free(c);
    for(i=0; i<=n_seq;i++){
      free(S1[i]);
      free(S2[i]);
    }
    free(S1); free(S2); //free(SS1); free(SS2);
    free(type);free(type2);
    return mfe;
    } else{
     struc = alibacktrack_XS(k_min, l_min,(const short int**)S1,(const short int**)S2,access_s1, access_s2);
    }
  int dx_5, dx_3, dy_5, dy_3,dGx,dGy,bonus_x, bonus_y,temp_dangle;
  dx_5=0; dx_3=0; dy_5=0; dy_3=0;dGx=0;dGy=0;bonus_x=0; bonus_y=0;temp_dangle=0;
  dGx =n_seq*(access_s1[i-k_min+1][i_pos]);dx_3=0; dx_5=0;bonus_x=0; 
  dGy = n_seq*access_s2[l_min-1][j_pos + (l_min-1)-1]; 
  mfe.tb=i_pos -9 - i + k_min -1 -dx_5;
  mfe.te=i_pos -9 -1 + dx_3;
  mfe.qb=j_pos -9 -1 - dy_5;
  mfe.qe=j_pos + l_min -3 -9 + dy_3;
  mfe.ddG=(double) (Emin *0.01);
  mfe.dG1=(double) (dGx*(0.01));
  mfe.dG2=(double) (dGy*(0.01));
  mfe.energy= mfe.ddG - mfe.dG1 - mfe.dG2;
  mfe.structure = struc;
  for (i=0; i<=n3; i++) free(c[i]);
  free(c);
  for(i=0; i<=n_seq;i++){
    free(S1[i]);
    free(S2[i]);
  }
  free(S1); free(S2); free(type);free(type2);
  return mfe;
}

PRIVATE char *alibacktrack_XS(int i, int j, const short *S1[], const short *S2[], 
                              const int** access_s1, const int ** access_s2) {
  int n3,n4,k, l, *type, type2, E, traced, i0, j0,s,n_seq,psc;
  char *st1=NULL, *st2=NULL, *struc=NULL;
  
  n3 = (int) S1[0][0];
  n4 = (int) S2[0][0];
  for (s=0; S1[s]!=NULL; s++);
  n_seq = s;
  for (s=0; S2[s]!=NULL; s++);
  if (n_seq != s) nrerror("unequal number of sequences in alibacktrack()\n");

  st1 = (char *) space(sizeof(char)*(n3+1));
  st2 = (char *) space(sizeof(char)*(n4+1));
  type = (int *) space(n_seq*sizeof(int));
  
  i0=i;/*MAX2(i-1,1);*/j0=j;/*MIN2(j+1,n4);*/
  while (i<=n3-1 && j>=2) {
    E = c[i][j]; traced=0;
    st1[i-1] = '(';
    st2[j-1] = ')'; 
    for (s=0; s<n_seq; s++) {
      type[s] = pair[S1[s][i]][S2[s][j]];
    }
    psc = covscore(type,n_seq);
    for (s=0; s<n_seq; s++) if (type[s]==0) type[s] = 7;
    E += psc;
    for (k=i+1; k<=n3 && k>i-MAXLOOP-2; k++) {
      for (l=j-1; l>=1; l--) {
        int LE;
        if (i-k+l-j-2>MAXLOOP) break;
        for (s=LE=0; s<n_seq; s++) {
          type2 = pair[S1[s][k]][S2[s][l]];
          if (type2==0) type2=7;
          LE += E_IntLoop(k-i-1, j-l-1, type[s], rtype[type2], 
                          S1[s][i+1], S2[s][j-1], S1[s][k-1], S2[s][l+1],P);
        }
        if (E == c[k][l]+LE) {
          traced=1; 
          i=k; j=l;
          break;
        }
      }
      if (traced) break;
    }
    if (!traced) { 
      for (s=0; s<n_seq; s++) {
        if (type[s]>2) E -= P->TerminalAU;
      }
      break;
      if (E != n_seq*P->DuplexInit) {
        nrerror("backtrack failed in fold duplex bal");
      } else break;
    }
  }
  struc = (char *) space(i-i0+1+j0-j+1+2);
  for (k=MAX2(i0,1); k<=i; k++) if (!st1[k-1]) st1[k-1] = '.';
  for (k=j; k<=j0; k++) if (!st2[k-1]) st2[k-1] = '.';
  strcpy(struc, st1+MAX2(i0-1,0)); strcat(struc, "&"); 
  strcat(struc, st2+j-1);
  free(st1); 
  free(st2);
  free(type);
  return struc;
}

duplexT** aliLduplexfold_XS(const char*s1[], const char* s2[], const int **access_s1, const int **access_s2, const int threshold, const int alignment_length, const int delta,  const int fast,const int il_a, const int il_b, const int b_a, const int b_b)
{
  short **S1, **S2;
  int *type,type2;
  int i, j,s,n_seq;
  s=0;
  int bopen=b_b;
  int bext=b_a;
  int iopen=il_b;
  int iext_s=2*(il_a);//iext_s 2 nt nucleotide extension of interior loop, on i and j side
  int iext_ass=50+il_a;//iext_ass assymetric extension of interior loop, either on i or on j side.
  int min_colonne=INF; //enthaelt das maximum einer kolonne
  int i_length;
  int max_pos;//get position of the best hit
  int max_pos_j;
  int temp;
  int min_j_colonne;
  int max=INF;
  //FOLLOWING NEXT 4 LINE DEFINES AN ARRAY CONTAINING POSITION OF THE SUBOPT IN S1
  int *position; //contains the position of the hits with energy > E
  int *position_j;
  int maxPenalty[4];
  
  n1 = (int) strlen(s1[0]);
  n2 = (int) strlen(s2[0]);
  for (s=0; s1[s]; s++);
  n_seq = s;
  for (s=0; s2[s]; s++);
  if (n_seq != s) nrerror("unequal number of sequences in aliduplexfold()\n");

  position = (int *) space ((delta+(n1)+4+delta) * sizeof(int));
  position_j= (int *) space ((delta+(n1)+4+delta) * sizeof(int));
  
  /**
  *** extension penalty, computed only once, further reduce the computation time
  **/
  if ((!P) || (fabs(P->temperature - temperature)>1e-6)){
    update_dfold_params();
  }
  maxPenalty[0]=(int) -1*P->stack[2][2]/2;
  maxPenalty[1]=(int) -1*P->stack[2][2];
  maxPenalty[2]=(int) -3*P->stack[2][2]/2;
  maxPenalty[3]=(int) -2*P->stack[2][2];
   

  lc   = (int**) space (sizeof(int *) * 5);
  lin  = (int**) space (sizeof(int *) * 5);
  lbx  = (int**) space (sizeof(int *) * 5);
  lby  = (int**) space (sizeof(int *) * 5);
  linx = (int**) space (sizeof(int *) * 5);
  liny = (int**) space (sizeof(int *) * 5);

  for (i=0; i<=4; i++){
    lc[i]  = (int *) space(sizeof(int) * (n2+5));  
    lin[i] = (int *) space(sizeof(int) * (n2+5));  
    lbx[i] = (int *) space(sizeof(int) * (n2+5));  
    lby[i] = (int *) space(sizeof(int) * (n2+5));  
    linx[i]= (int *) space(sizeof(int) * (n2+5));  
    liny[i]= (int *) space(sizeof(int) * (n2+5));  
  }

  
  S1 = (short **) space((n_seq+1)*sizeof(short *));
  S2 = (short **) space((n_seq+1)*sizeof(short *));
  for (s=0; s<n_seq; s++) {
    if (strlen(s1[s]) != n1) nrerror("uneqal seqence lengths");
    if (strlen(s2[s]) != n2) nrerror("uneqal seqence lengths");
    S1[s] = encode_seq(s1[s]);
    S2[s] = encode_seq(s2[s]);    
  }
  type = (int *) space(n_seq*sizeof(int));
  /**
  *** array initialization
  **/
    
  for(j=n2+4;j>=0;j--) {
    lbx[0][j]=lbx[1][j]=lbx[2][j]=lbx[3][j]    = lbx[4][j] =INF;
    lin[0][j]=lin[1][j]=lin[2][j]=lin[3][j]    = lin[4][j] =INF;
    lc[0][j] =lc[1][j] =lc[2][j] = lc[3][j]    =  lc[4][j] =INF;
    lby[0][j]=lby[1][j]=lby[2][j]=lby[3][j]    = lby[4][j] =INF;
    liny[0][j]=liny[1][j]=liny[2][j]=liny[3][j]=liny[4][j]=INF;
    linx[0][j]=linx[1][j]=linx[2][j]=linx[3][j]=linx[4][j]=INF;
  }
  i=10;
  i_length= n1 - 9 ;
  int di1,di2,di3,di4; //contains accessibility penalty
  while(i < i_length) {
    int idx=i%5;
    int idx_1=(i-1)%5;
    int idx_2=(i-2)%5;
    int idx_3=(i-3)%5;
    int idx_4=(i-4)%5;
    
    di1 = access_s1[5][i]   - access_s1[4][i-1];           
    di2 = access_s1[5][i-1] - access_s1[4][i-2] + di1;
    di3 = access_s1[5][i-2] - access_s1[4][i-3] + di2;
    di4 = access_s1[5][i-3] - access_s1[4][i-4] + di3;
    di1=MIN2(di1,maxPenalty[0]);
    di2=MIN2(di2,maxPenalty[1]);
    di3=MIN2(di3,maxPenalty[2]);
    di4=MIN2(di3,maxPenalty[3]);
    j=n2-9;
    while (9 < --j) {
      int dj1,dj2,dj3,dj4;
      dj1 = access_s2[5][j+4] - access_s2[4][j+4];        
      dj2 = access_s2[5][j+5] - access_s2[4][j+5] + dj1;
      dj3 = access_s2[5][j+6] - access_s2[4][j+6] + dj2;
      dj4 = access_s2[5][j+7] - access_s2[4][j+7] + dj3;
      dj1=MIN2(dj1,maxPenalty[0]);
      dj2=MIN2(dj2,maxPenalty[1]);
      dj3=MIN2(dj3,maxPenalty[2]);
      dj4=MIN2(dj4,maxPenalty[3]);
      int psc;
      for (s=0; s<n_seq; s++) { //initialize type1
        type[s] = pair[S1[s][i]][S2[s][j]];
      }
      psc = covscore(type, n_seq); 
      for (s=0; s<n_seq; s++) if (type[s]==0) type[s]=7;
      lc[idx][j] = ((psc >= MINPSCORE) ? n_seq*(P->DuplexInit + access_s1[1][i] + access_s2[1][j]) : INF);
      int c_stack, c_10, c_01, c_11, c_22, c_21, c_12, c_23, c_32, c_in, c_in2x, c_in2y, c_bx, c_by, c_inx, c_iny;  //matrix c
      int in,  in_x, in_y, in_xy; // in begin, in_x assymetric, in_y assymetric, in_xy symetric;
      int inx, inx_x;
      int iny, iny_y;
      int bx,  bx_x; 
      int by,  by_y; 
      in=lc[idx_1][j+1]; in_x=lin[idx_1][j]; in_y=lin[idx][j+1]; in_xy=lin[idx_1][j+1];
      inx=lc[idx_1][j+1]; inx_x=linx[idx_1][j];
      iny=lc[idx_1][j+1]; iny_y=liny[idx][j+1];
      bx=lc[idx_1][j]; bx_x=lbx[idx_1][j]; 
      by=lc[idx][j+1]; by_y=lby[idx][j+1];
      c_stack=lc[idx_1][j+1]; c_01=lc[idx_1][j+2];c_10=lc[idx_2][j+1]; 
      c_12=lc[idx_2][j+3];c_21=lc[idx_3][j+2];c_11=lc[idx_2][j+2];
      c_22=lc[idx_3][j+3];c_32=lc[idx_4][j+3];c_23=lc[idx_3][j+4];
      c_in=lin[idx_3][j+3];c_in2x=lin[idx_4][j+2];c_in2y=lin[idx_2][j+4];
      c_inx=linx[idx_3][j+1]; c_iny=liny[idx_1][j+3];
      c_bx=lbx[idx_2][j+1];c_by=lby[idx_1][j+2];
      for (s=0; s<n_seq; s++) {
        type2 = pair[S2[s][j+1]][S1[s][i-1]];
        in   +=P->mismatchI[type2][S2[s][j]][S1[s][i]]+iopen+iext_s+di1+dj1;
        in_x +=iext_ass+di1;
        in_y +=iext_ass+dj1;
        in_xy+=iext_s+di1+dj1;
        inx  +=P->mismatch1nI[type2][S2[s][j]][S1[s][i]]+iopen+iext_s+di1+dj1;
        inx_x+=iext_ass+di1;
        iny  +=P->mismatch1nI[type2][S2[s][j]][S1[s][i]]+iopen+iext_s+di1+dj1;
        iny_y+=iext_ass+dj1;
        type2=pair[S2[s][j]][S1[s][i-1]];   
        bx   +=bopen+bext+(type2>2?P->TerminalAU:0)+di1;
        bx_x +=bext+di1;
        type2=pair[S2[s][j+1]][S1[s][i]];
        by   +=bopen+bext+(type2>2?P->TerminalAU:0)+dj1;
        by_y +=bext+dj1;
      }
      lin[idx][j]=MIN2(in, MIN2(in_x, MIN2(in_y, in_xy)));       
      linx[idx][j]=MIN2(inx_x, inx);
      liny[idx][j]=MIN2(iny_y, iny);
      lby[idx][j]=MIN2(by, by_y);
      lbx[idx][j]=MIN2(bx, bx_x);
      if (psc<MINPSCORE) continue;  
      for (s=0; s<n_seq; s++) {
        lc[idx][j]+=E_ExtLoop(type[s], S1[s][i-1],S2[s][j+1],P);
      }
      for (s=0; s<n_seq; s++) {
        type2=pair[S1[s][i-1]][S2[s][j+1]];if (type2==0) type2=7;
        c_stack+=E_IntLoop(0,0,type2, rtype[type[s]],S1[s][i], S2[s][j], S1[s][i-1], S2[s][j+1], P)+di1+dj1;
        type2=pair[S1[s][i-1]][S2[s][j+2]];if (type2==0) type2=7;
        c_01   +=E_IntLoop(0,1,type2, rtype[type[s]],S1[s][i], S2[s][j+1], S1[s][i-1], S2[s][j+1], P)+di1+dj2;
        type2=pair[S1[s][i-2]][S2[s][j+1]];if (type2==0) type2=7;
        c_10   +=E_IntLoop(1,0,type2, rtype[type[s]],S1[s][i-1], S2[s][j], S1[s][i-1], S2[s][j+1], P)+di2+dj1;
        type2=pair[S1[s][i-2]][S2[s][j+2]];if (type2==0) type2=7;
        c_11   +=E_IntLoop(1,1,type2, rtype[type[s]],S1[s][i-1], S2[s][j+1], S1[s][i-1], S2[s][j+1], P)+di2+dj2;
        type2=pair[S1[s][i-3]][S2[s][j+3]];if (type2==0) type2=7;
        c_22   +=E_IntLoop(2,2,type2, rtype[type[s]],S1[s][i-2], S2[s][j+2], S1[s][i-1], S2[s][j+1], P)+di3+dj3;
        type2= pair[S1[s][i-3]][S2[s][j+2]];if (type2==0) type2=7;
        c_21   +=E_IntLoop(2,1,type2, rtype[type[s]],S1[s][i-2], S2[s][j+1], S1[s][i-1], S2[s][j+1], P)+di3+dj2;
        type2= pair[S1[s][i-2]][S2[s][j+3]];if (type2==0) type2=7;
        c_12   +=E_IntLoop(1,2,type2, rtype[type[s]],S1[s][i-1], S2[s][j+2], S1[s][i-1], S2[s][j+1], P)+di2+dj3;
        type2 = pair[S1[s][i-4]][S2[s][j+3]];if (type2==0) type2=7;
        c_32   +=E_IntLoop(3,2,type2, rtype[type[s]],S1[s][i-3], S2[s][j+2], S1[s][i-1], S2[s][j+1], P)+di4+dj3;
        type2 = pair[S1[s][i-3]][S2[s][j+4]];if (type2==0) type2=7;
        c_23   +=E_IntLoop(2,3,type2, rtype[type[s]],S1[s][i-2], S2[s][j+3], S1[s][i-1], S2[s][j+1], P)+dj4+di3;
        c_in   +=P->mismatchI[rtype[type[s]]][S1[s][i-1]][S2[s][j+1]]+di3+dj3+2*iext_s;
        c_in2x +=P->mismatchI[rtype[type[s]]][S1[s][i-1]][S2[s][j+1]]+iext_s+2*iext_ass+di4+dj2;
        c_in2y +=P->mismatchI[rtype[type[s]]][S1[s][i-1]][S2[s][j+1]]+iext_s+2*iext_ass+di2+dj4;
        c_inx  +=P->mismatch1nI[rtype[type[s]]][S1[s][i-1]][S2[s][j+1]]+iext_ass+iext_ass+di3+dj1;
        c_iny  +=P->mismatch1nI[rtype[type[s]]][S1[s][i-1]][S2[s][j+1]]+iext_ass+iext_ass+di1+dj3;
        int bAU;
        bAU=(type[s]>2?P->TerminalAU:0);
        c_bx   +=di2+dj1+bext+bAU;
        c_by   +=di1+dj2+bext+bAU;
      }
      lc[idx][j] =MIN2(lc[idx][j], 
                       MIN2(c_stack, 
                            MIN2(c_10, 
                                 MIN2(c_01, 
                                      MIN2(c_11, 
                                           MIN2(c_21, 
                                                MIN2(c_12, 
                                                     MIN2(c_22, 
                                                          MIN2(c_23, 
                                                               MIN2(c_32, 
                                                                    MIN2(c_bx, 
                                                                         MIN2(c_by, 
                                                                              MIN2(c_in,
                                                                                   MIN2(c_in2x,
                                                                                        MIN2(c_in2y,
                                                                                             MIN2(c_inx,c_iny)
                                                                                             )
                                                                                        )
                                                                                   )
                                                                              )
                                                                         )
                                                                    )
                                                               )
                                                          )
                                                     )
                                                )
                                           )
                                      )
                                 )
                            )
                       );
      lc[idx][j]-=psc;
      temp=lc[idx][j];

      for (s=0; s<n_seq; s++) {
        temp+=E_ExtLoop(rtype[type[s]], S2[s][j-1],S1[s][i+1],P);
      }
      if(min_colonne > temp){
        min_colonne=temp;
        min_j_colonne=j;
      }
    }
    if(max>=min_colonne){
            max=min_colonne;
            max_pos=i;
        max_pos_j=min_j_colonne;
    }
    position[i+delta]=min_colonne;min_colonne=INF;
    position_j[i+delta]=min_j_colonne;
    i++;
  }
  //printf("MAX: %d",max);
  for (s=0; s<n_seq; s++) {free(S1[s]);free(S2[s]);}
  free(S1); free(S2); 
  if(max<threshold){
    alifind_max_XS(position, position_j, delta, threshold, alignment_length, s1, s2,access_s1,access_s2, fast);
  }
  aliplot_max_XS(max, max_pos, max_pos_j,alignment_length, s1, s2, access_s1, access_s2, fast); 
  for (i=0; i<=4; i++) {free(lc[i]);free(lin[i]);free(lbx[i]);free(lby[i]);free(linx[i]);free(liny[i]);}
  //free(lc[0]);free(lin[0]);free(lbx[0]);free(lby[0]);free(linx[0]);free(liny[0]);
  free(lc);free(lin);free(lbx);free(lby);free(linx);free(liny);
  free(position);
  free(position_j);
  free(type);
  return NULL;
}




PRIVATE void alifind_max_XS(const int *position, const int *position_j,
                                const int delta, const int threshold, const int alignment_length, 
                                const char *s1[], const char *s2[], 
                                const int **access_s1, const int **access_s2, const int fast){

  int n_seq=0;
  for (n_seq=0; s1[n_seq]!=NULL; n_seq++);
  int pos=n1-9;
  if(fast==1){
    while(10 < pos--){
      int temp_min=0;
      if(position[pos+delta]<(threshold)){                                        
        int search_range;                                                   
        search_range=delta+1;                                               
        while(--search_range){   
          if(position[pos+delta-search_range]<=position[pos+delta-temp_min]){
            temp_min=search_range;                                          
          }                                                                 
        }                                                                   
        pos-=temp_min;                                                      
        int max_pos_j;                                                      
        max_pos_j=position_j[pos+delta];
        int max;
        max=position[pos+delta];
        //        printf("target upper bound %d: query lower bound %d  (%5.2f) \n", pos-10, max_pos_j-10, ((double)max)/100);
        pos=MAX2(10,pos-delta);
      }
    }
  }
  else{
    pos=n1-9;
    while( pos-- > 10 ) {
      //printf("delta %d position:%d value:%d\n", delta, pos, position[pos]);
      int temp_min=0;
      if(position[pos+delta]<(threshold)){
        int search_range;
        search_range=delta+1;
        while(--search_range){
           if(position[pos+delta-search_range]<=position[pos+delta-temp_min]){
            temp_min=search_range;
          }
        }
        pos-=temp_min; 
        int max_pos_j;
        max_pos_j=position_j[pos+delta]; //position on j
        //printf("%d %d %d\n", pos, max_pos_j,position[pos+delta]);
        int begin_t=MAX2(9, pos-alignment_length); //only get the position that binds..
        int end_t  =MIN2(n1-10, pos);              //..no dangles
        int begin_q=MAX2(9, max_pos_j-2);
        int end_q  =MIN2(n2-10, max_pos_j+alignment_length-2);
        char **s3, **s4;
        s3 = (char**) space(sizeof(char*)*(n_seq+1));
        s4 = (char**) space(sizeof(char*)*(n_seq+1));
        int i;
        for(i=0; i<n_seq; i++){
          s3[i] = (char*) space(sizeof(char)*(end_t-begin_t+2));
          s4[i] = (char*) space(sizeof(char)*(end_q-begin_q+2));
          strncpy(s3[i], (s1[i]+begin_t), end_t - begin_t +1);
          strncpy(s4[i], (s2[i]+begin_q), end_q - begin_q +1);
          s3[i][end_t - begin_t +1]='\0';
          s4[i][end_q - begin_q +1]='\0';
        }
        duplexT test;
        test = aliduplexfold_XS((const char**) s3, (const char**) s4, access_s1, access_s2,pos, max_pos_j,threshold);
        //        printf("position %d approximation %d test %d threshold %d\n", pos, position[pos+delta], (int)test.energy,(int)(threshold/n_seq));
        if(test.energy*100  < (int) (threshold/n_seq)){
        printf("%s %3d,%-3d: %3d,%-3d (%5.2f = %5.2f + %5.2f + %5.2f)\n", test.structure,
               test.tb,test.te,test.qb,test.qe, test.ddG/n_seq, test.energy/n_seq, test.dG1/n_seq, test.dG2/n_seq);
        free(test.structure);
        pos=MAX2(10,pos-delta);
        }
        for(i=0;i<n_seq;i++){
          free(s3[i]);free(s4[i]);
        }
        free(s3);free(s4);
        //free(test.structure);
      }
    }
  }
}

PRIVATE void aliplot_max_XS(const int max, const int max_pos, const int max_pos_j,const int alignment_length, 
                            const char *s1[], const char* s2[],const int **access_s1, const int **access_s2, 
                            const int fast){

  int n_seq;
  for (n_seq=0; !(s1[n_seq]==NULL); n_seq++);
  n1 = strlen(s1[0]); //get length of alignment
  n2 = strlen(s2[0]); //get length of alignme
  
  if(fast){
    printf("target upper bound %d: query lower bound %d (%5.2f)\n", 
           max_pos-10, max_pos_j-10, (double) ((double)max)/(100*n_seq));
  }
  else{
    int begin_t=MAX2(9, max_pos-alignment_length); //only get the position that binds..
    int end_t  =MIN2(n1-10, max_pos);              //..no dangles
    int begin_q=MAX2(9, max_pos_j-2);
    int end_q  =MIN2(n2-10, max_pos_j+alignment_length-2);
    char **s3, **s4;
    s3 = (char**) space(sizeof(char*)*(n_seq+1));
    s4 = (char**) space(sizeof(char*)*(n_seq+1));
    int i;
    for(i=0; i<n_seq; i++){
      s3[i] = (char*) space(sizeof(char)*(end_t-begin_t+2));
      s4[i] = (char*) space(sizeof(char)*(end_q-begin_q+2));
      strncpy(s3[i], (s1[i]+begin_t), end_t - begin_t +1);
      strncpy(s4[i], (s2[i]+begin_q), end_q - begin_q +1);
      s3[i][end_t - begin_t +1]='\0';
      s4[i][end_q - begin_q +1]='\0';
    }
    duplexT test;
    test = aliduplexfold_XS((const char**) s3, (const char**) s4, access_s1, access_s2,max_pos, max_pos_j,INF);
    printf("%s %3d,%-3d: %3d,%-3d (%5.2f = %5.2f + %5.2f + %5.2f)\n", test.structure,
           test.tb,test.te,test.qb,test.qe, test.ddG/n_seq, test.energy/n_seq, test.dG1/n_seq, test.dG2/n_seq);
    free(test.structure);
    for(i=0;i<n_seq;i++){
      free(s3[i]);free(s4[i]);
    }
    free(s3);free(s4);
  }
}

PRIVATE int covscore(const int *types, int n_seq) {
  /* calculate co-variance bonus for a pair depending on  */
  /* compensatory/consistent mutations and incompatible seqs */
  /* should be 0 for conserved pairs, >0 for good pairs      */
#define NONE -10000 /* score for forbidden pairs */
  int k,l,s,score, pscore;
  int dm[7][7]={{0,0,0,0,0,0,0}, /* hamming distance between pairs */
                {0,0,2,2,1,2,2} /* CG */,
                {0,2,0,1,2,2,2} /* GC */,
                {0,2,1,0,2,1,2} /* GU */,
                {0,1,2,2,0,2,1} /* UG */,
                {0,2,2,1,2,0,2} /* AU */,
                {0,2,2,2,1,2,0} /* UA */};
  
  int pfreq[8]={0,0,0,0,0,0,0,0};
  for (s=0; s<n_seq; s++)
    pfreq[types[s]]++;

  if (pfreq[0]*2>n_seq) return NONE;
  for (k=1,score=0; k<=6; k++) /* ignore pairtype 7 (gap-gap) */
    for (l=k+1; l<=6; l++)
      /* scores for replacements between pairtypes    */
      /* consistent or compensatory mutations score 1 or 2  */
      score += pfreq[k]*pfreq[l]*dm[k][l];
  
  /* counter examples score -1, gap-gap scores -0.25   */
  pscore = cv_fact * 
    ((UNIT*score)/n_seq - nc_fact*UNIT*(pfreq[0] + pfreq[7]*0.25));
  return pscore;
}


PRIVATE void update_dfold_params(void)
{
  if(P) free(P);
  P = scale_parameters();
  make_pair_matrix();
}



PRIVATE short * encode_seq(const char *sequence) {
  unsigned int i,l;
  short *S;
  l = strlen(sequence);
  S = (short *) space(sizeof(short)*(l+2));
  S[0] = (short) l;

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

  /* for circular folding add first base at position n+1 */
  S[l+1] = S[1];

  return S;
}