#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <stdarg.h>
#include <string.h>

#include "io_lib_header.h"
#include "util_lib_header.h"
#include "define_header.h"
#include "dp_lib_header.h"


int cl2pair_list_ecl ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in);



/*******************************************************************************/
/*                idscore_pairseq: measure the % id without delivering thze aln*/                                                   
/*                                                                             */
/*	makes DP between the the ns[0] sequences and the ns[1]                 */
/*                                                                             */
/*	for MODE, see the function  get_dp_cost                                */
/*******************************************************************************/
int idscore_pairseq (char *s1, char *s2, int gop, int gep, int **m, char *comp_mode)
{
  int **I, **D, **M, *P;
  int i, j,l1, l2, l,score, id, igop,match;


  l1=strlen (s1); l2=strlen (s2);
  lower_string (s1); lower_string (s2);
  
  I=declare_int (6,l2+1);D=declare_int (6,l2+1);M=declare_int (6,l2+1);
  for (j=0; j<=l2; j++)
    {
      D[0][j]=gep*j;M[0][j]=2*gep*j;D[4][j]=0;
    }
  
  for (i=1; i<=l1; i++)
    {

      I[1][0]=i*gep;
      M[1][0]=2*i*gep;
      
      for (j=1; j<=l2; j++)
	{
	  score=m[s1[i-1]-'a'][s2[j-1]-'a'];	  
	  id=(s1[i-1]==s2[j-1])?1:0;
	  
	  igop=(i==l1 || j==l2)?0:gop;

	  if   ((D[0][j]+gep)>(M[0][j]+igop+gep))   {D[1][j]=D[0][j]+gep;      D[3][j]=D[2][j];        D[5][j]=D[4][j];}
	  else                                      {D[1][j]=M[0][j]+igop+gep; D[3][j]=M[2][j];        D[5][j]=M[4][j];}
	  
	  if ( (I[1][j-1]+gep)>(M[1][j-1]+igop+gep)){I[1][j]=I[1][j-1]+gep;      I[3][j]=I[3][j-1];    I[5][j]=I[5][j-1];}
	  else                                      {I[1][j]=M[1][j-1]+igop+gep; I[3][j]=M[3][j-1];    I[5][j]=M[5][j-1];}
	  
	  match=M[0][j-1]+score;
	  if (I[1][j]>match && I[1][j]>D[1][j])     {M[1][j]=I[1][j]           ; M[3][j]=I[3][j];      M[5][j]=I[5][j];}
	  else if (D[1][j]>match)                   {M[1][j]=D[1][j]           ; M[3][j]=D[3][j];      M[5][j]=D[5][j];}
	  else                                      {M[1][j]=match             ; M[3][j]=M[2][j-1]+id; M[5][j]=M[4][j-1]+1;}
	}
      P=I[0]; I[0]=I[1]; I[1]=P;
      P=I[2]; I[2]=I[3]; I[3]=P;
      P=I[4]; I[4]=I[5]; I[5]=P;
      
      P=D[0]; D[0]=D[1]; D[1]=P;
      P=D[2]; D[2]=D[3]; D[3]=P;
      P=D[4]; D[4]=D[5]; D[5]=P;
      
      P=M[0]; M[0]=M[1]; M[1]=P;
      P=M[2]; M[2]=M[3]; M[3]=P;
      P=M[4]; M[4]=M[5]; M[5]=P;
    }
 

  

  if ( strstr (comp_mode, "sim2"))
    {
       l=MIN(l1,l2);
       score=(l==0)?0:(M[2][l2]*100)/l;
    }
  else if ( strstr (comp_mode, "sim3"))
    {
       l=MAX(l1,l2);
       score=(l==0)?0:(M[2][l2]*100)/l;
    }
  else if ( strstr (comp_mode, "cov"))
    {
      l=MAX(l1,l2);
      score=(l==0)?0:((M[4][l2]*100)/l);
    }
  else
    {
      //default: simple sim
      l=M[4][l2];
      score=(l==0)?0:(M[2][l2]*100)/l;
    }      
  
  free_int (I, -1);
  free_int (D, -1);
  free_int (M, -1);
  
  return score;
}
	      
int test_pair_wise (Alignment *A, int *ns, int **l_s, Constraint_list *CL)
{
  int a,l0, l1, n;
  char buf[VERY_LONG_STRING];
  char *gap, *seq;
  
  l0=strlen (A->seq_al[l_s[0][0]]);
  l1=strlen (A->seq_al[l_s[1][0]]);
  
  n=(l0<5)?l0/2:5;
  gap=generate_null(l1-n);
  for (a=0;a<ns[0]; a++)
    {
      seq=A->seq_al[l_s[0][a]];
      sprintf (buf, "%s%s",seq, gap);
      sprintf (seq, "%s", buf);
    }
  vfree (gap);
  gap=generate_null(l0-n);
  
  for (a=0;a<ns[1]; a++)
    {
      seq=A->seq_al[l_s[1][a]];
      sprintf (buf, "%s%s",seq, gap);
      sprintf (seq, "%s", buf);
    }
  vfree(gap);
  

  A->len_aln=strlen (A->seq_al[l_s[0][0]]);
  A->score=A->score_aln=100;
  return 100;
}  

int idscore_pair_wise (Alignment *A,int*ns, int **l_s,Constraint_list *CL)
{
  
  A->score_aln=A->score=idscore_pairseq (A->seq_al[l_s[0][0]], A->seq_al[l_s[1][0]], CL->gop, CL->gep,CL->M, "sim3");
  return A->score_aln;
}
int dp_max (int *trace, int n, ...);
int dp_max (int *trace, int n, ...)
{
  va_list ap;
  int a, v, t, best_v=0;
  
  va_start (ap, n);
  for (a=0; a< n; a++)
    {
      t=va_arg (ap, int);
      v=va_arg (ap, int);

      if (a==0)
	{
	  best_v=v;
	  trace[0]=t;
	}
      else
	{
	  if (v>best_v)
	    {
	      best_v=v;
	      trace[0]=t;
	    }
	}
    }
 
  return best_v;
}
int is_tied (int *trace, int n, ...);
int is_tied(int *trace, int n, ...)
{
  va_list ap;
  int a, v, t, best_v=0;
  int nties=0;
  
  va_start (ap, n);
  for (a=0; a< n; a++)
    {
      t=va_arg (ap, int);
      v=va_arg (ap, int);

      if (a==0)
	{
	  best_v=v;
	  trace[0]=t;
	}
      else
	{
	  if (v>best_v)
	    {
	      best_v=v;
	      trace[0]=t;
	    }
	}
    }
  va_end(ap);
  va_start (ap,n);
  for (a=0; a<n; a++)
    {
      t=va_arg (ap, int);
      v=va_arg (ap, int);
      if (v==best_v && trace[0]!=t)
	nties++;
    }
  va_end (ap);
  return nties;
}

void display_mat (int **M, int l1, int l2, char *title);
void display_mat (int **M, int l1, int l2, char *title)
{
  int a, b;
  
  fprintf ( stdout, "\n\nTitle %s\n", title);
  for ( a=0; a<=l1; a++)
    {
      fprintf ( stdout, "\n");
      for ( b=0; b<=l2; b++)
	fprintf ( stdout, "%3d ", M[a][b]);
    }
}
int glocal_pair_wise ( Alignment *A, int *ns, int **l_s, Constraint_list *CL)
{
  int ***t, ***m;
  int i,j, l1, l2, n, sub, trace,ntrace, a, b, c, score;
  int gop,rgop,tgop, gep, unmatch;
  int M1, M2, I1, D1, LEN;
  char **al, *char_buf, **aln;
  int **pos0;
  

  l1=strlen (A->seq_al[l_s[0][0]]);
  l2=strlen (A->seq_al[l_s[1][0]]);

  n=1;
  M1=n++;D1=n++;I1=n++;M2=n++;
  t=declare_arrayN(3, sizeof (int),n, l1+1, l2+1);
  m=declare_arrayN(3, sizeof (int),n, l1+1, l2+1);
  
  
  gop=CL->gop*SCORE_K;
  gep=CL->gep*SCORE_K;
  tgop=gop;
  unmatch=gep;
  
  pos0=aln2pos_simple ( A,-1, ns, l_s);
 
  
  for (j=1; j<=l2; j++)
    {
      m[D1][0][j]=gep*j;

      m[M1][0][j]=2*gep*j;
      m[M2][0][j]=4*gep*j;
    }
  
  
  for (i=1; i<=l1; i++)
    {
      m[I1][i][0]=i*gep;
      m[M2][i][0]=4*i*gep;
      m[M1][i][0]=2*i*gep;
                 
      for ( j=1; j<=l2; j++)
	{
	  rgop=(i==l1 || j==1)?0:gop;
	  rgop=gop;
	  sub=(CL->get_dp_cost) (A, pos0, ns[0], l_s[0], i-1, pos0, ns[1], l_s[1],j-1,CL);	
	  m[M1][i][j]=dp_max (&trace,4,M1,m[M1][i-1][j-1],I1, m[I1][i-1][j-1],D1,m[D1][i-1][j-1],M2,m[M2][i-1][j-1])+sub;
	  t[M1][i][j]=trace;
	  
	  m[D1][i][j]=dp_max (&trace,3, M1,m[M1][i][j-1]+rgop,D1, m[D1][i][j-1]+gep, M2, m[M2][i][j-1]);
	  t[D1][i][j]=trace;
	  
	  m[I1][i][j]=dp_max (&trace,3, M1,m[M1][i-1][j]+rgop, I1, m[I1][i-1][j]+gep, M2, m[M2][i-1][j]);
	  t[I1][i][j]=trace;

	  m[M2][i][j]=dp_max (&trace,4,M1,m[M1][i-1][j-1]+tgop,I1, m[I1][i-1][j-1]+tgop,D1,m[D1][i-1][j-1]+tgop,M2,m[M2][i-1][j-1])+unmatch;
	  t[M2][i][j]=trace;
	  
	}
	  
    }
  score=dp_max (&trace,4, M1,m[M1][l1][l2],D1,m[D1][l1][l2],I1, m[I1][l1][l2],M2,m[M2][l1][l2]);
  LEN=0;i=l1;j=l2;
  al=declare_char (2, l1+l2+1);
  

  trace=t[trace][i][j];
  while (!(i==0 &&j==0))
    {
  
      ntrace=t[trace][i][j];
      if (i==0)
	{
	  al[0][LEN]=0;
	  al[1][LEN]=1;
	  j--;
	  LEN++;
	}
      else if ( j==0)
	{
	  al[0][LEN]=1;
	  al[1][LEN]=0;
	  i--;
	  LEN++;
	}
      else if ( trace==M1)
	{
	  al[0][LEN]=1;
	  al[1][LEN]=1;
	  i--; j--;
	  LEN++;
	}
      else if ( trace==M2)
	{
	  al[0][LEN]=1;
	  al[1][LEN]=0;
	  LEN++;

	  al[0][LEN]=0;
	  al[1][LEN]=1;
	  LEN++;

	  i--; j--;
	  
	}
      else if ( trace==D1)
	{
	  al[0][LEN]=0;
	  al[1][LEN]=1;
	  j--;
	  LEN++;
	}
      else if ( trace == I1)
	{
	  al[0][LEN]=1;
	  al[1][LEN]=0;
	  i--;
	  LEN++;
	}
      trace=ntrace;     
      
    }
  
  invert_list_char ( al[0], LEN);
  invert_list_char ( al[1], LEN);	
  if ( A->declared_len<=LEN)A=realloc_aln2  ( A,A->max_n_seq, 2*LEN);	
  
  aln=A->seq_al;
  char_buf= vcalloc (LEN+1, sizeof (char));	
  for ( c=0; c< 2; c++)
    {
      for ( a=0; a< ns[c]; a++) 
	{		
	  int ch=0;
	  for ( b=0; b< LEN; b++)
	    {		   
	      if (al[c][b]==1)
		char_buf[b]=aln[l_s[c][a]][ch++];
	      else
		char_buf[b]='-';
	    }
	  char_buf[b]='\0';
	  sprintf (aln[l_s[c][a]],"%s", char_buf);
	}
    }
  
  
  A->len_aln=LEN;
  A->nseq=ns[0]+ns[1];
  free_arrayN((void *)m, 3);
  free_arrayN((void *)t, 3);
  vfree (char_buf);
  free_char (al, -1);
  return score;
}

int ** aln2local_penalties (Alignment *A, int n, int *ls, Constraint_list *CL, int **lg);
int ** aln2local_penalties (Alignment *A, int n, int *ls, Constraint_list *CL, int **lg)
{
  //adapted from gap_count in MAFFT V 5.5
  int p,s,l, c1, c2;
  int gep,gop;
  int open=3, close=4, gap=5;
  
  gop=CL->gop*SCORE_K;
  gep=CL->gep*SCORE_K;
  
  l=strlen (A->seq_al[ls[0]]);
  
  if (!lg)
    {
      lg=declare_int (6, l);
    }
  
  if ( read_array_size_new (lg[0])<l)
    {
      free_int (lg, -1);
      lg=declare_int (6, l);
    }
  
  for( s=0; s<n; s++ ) 
	{
	  c1='x';
	  for (p=0; p<l; p++)
	    {
	      c2=A->seq_al[ls[s]][p];
	      
	      if (c1!='-' && c2=='-')lg[open][p]++;
	      if (c1=='-' && c2!='-')lg[close][p]++;
	      if ( c1=='-')lg[gap][p]++;
	      c1=c2;
	    }
	}
  
  for (p=0; p<l; p++)
    {
      float go, gc, nn;
      nn=n;
      go=lg[open ][p];
      gc=lg[close][p];
     
    
      lg[GOP][p]=0.5*(1-(go/nn))*gop;
      lg[GCP][p]=0.5*(1-(gc/nn))*gop;
      //Checked locacal gep => gives low quality results
      lg[GEP][p]=gep;//(1-((float)lg[gap][p]/(float)n))*gep;
      lg[open][p]=lg[close][p]=lg[gap][p]=0;
      
    }

  return lg;
}
int free_gotoh_pair_wise_lgp()
{
  return gotoh_pair_wise_lgp (NULL, NULL, NULL, NULL);
}
int gotoh_pair_wise_lgp ( Alignment *A, int *ns, int **l_s, Constraint_list *CL)
{
  int i,j, li, lj, n, sub, trace,ntrace, a, b, c, score;
  int I, J;
  int M1, I1, D1, LEN;
  char **al, *char_buf, **aln;
  int **pos0, **pos;
  Alignment *Aln;
  
  int gop[2], gcp[2], gep[2];
  static int ***gpl, ***t, ***m;
  static int max_li, max_lj;
  
  

  //gotoh_pair_wise ( A, ns, l_s,CL);
  //ungap_sub_aln (A, ns[0], l_s[0]);
  //ungap_sub_aln (A, ns[1], l_s[1]);
  
  if (!A)
    {
      free_arrayN((void**)gpl, 3);
      free_arrayN((void**)t, 3);
      free_arrayN((void**)m, 3);
      max_li=max_lj=0;
      return 0;
    }

  I=0;J=1;

  
  li=strlen (A->seq_al[l_s[I][0]]);
  lj=strlen (A->seq_al[l_s[J][0]]);
  
  if ( !gpl)gpl=vcalloc ( 2, sizeof (int**));
  gpl[I]=aln2local_penalties (A,ns[I], l_s[I], CL,gpl[I]);
  gpl[J]=aln2local_penalties (A,ns[J], l_s[J], CL,gpl[J]);
  
  
  n=1;
  M1=n++;D1=n++;I1=n++;
  
  if ( li>max_li ||lj>max_lj )
    {
      free_arrayN((void**)t, 3);
      free_arrayN((void**)m, 3);

     
      max_li=li;
      max_lj=lj;
      t=declare_arrayN(3, sizeof (int),n, max_li+1, max_lj+1);
      m=declare_arrayN(3, sizeof (int),n, max_li+1, max_lj+1);
      
    }
  pos0=aln2pos_simple ( A,-1, ns, l_s);
 
  //Compatibility with Macro
  Aln=A;
  pos=pos0;
  
  for (j=1; j<=lj; j++)
    {
      gep[J]=gpl[J][GEP][j-1];
      m[D1][0][j]=gep[J]*j;
      m[I1][0][j]=m[D1][0][j]-1;
      m[M1][0][j]=m[D1][0][j]-1;
    }
  
  //D1: gap in sequence I
  //I1: gap in sequence J
  
  
  for (i=1; i<=li; i++)
    {
      gep[I]=gpl[I][GEP][i-1];
      gop[I]=gpl[I][GOP][i-1];
      gcp[I]=gpl[I][GCP][i-1];
      
      m[I1][i][0]=i*gep[I];
      m[D1][i][0]= m[I1][i][0]-1;
      m[M1][i][0]= m[I1][i][0]-1;
                 
     
      
      gop[I]=(i==1 || i==li )?0:gop[I];
      gcp[I]=(i==1 || i==li )?0:gcp[I];
      
      
      for ( j=1; j<=lj; j++)
	{
	  
	  gep[J]=gpl[J][GEP][j-1];
	  gop[J]=gpl[J][GOP][j-1];
	  gcp[J]=gpl[J][GCP][j-1];
	  
	  //gep[J]=gep[I]=(gep[J]+gep[I])/2;
	  //gop[J]=gop[I]=(gop[J]+gop[I])/2;
	  //gcp[J]=gcp[I]=(gcp[J]+gcp[I])/2;
	  

	  gop[J]=(j==1 || j==lj )?0:gop[J];
	  gcp[J]=(j==1 || j==lj )?0:gcp[J];
	  
	  
	  //sub=(CL->get_dp_cost) (A, pos0, ns[0], l_s[0], i-1, pos0, ns[1], l_s[1],j-1,CL);	
	  sub=TC_SCORE((i-1), (j-1));
	  
	  m[M1][i][j]=dp_max (&trace,3,M1,m[M1][i-1][j-1],I1, m[I1][i-1][j-1]+gcp[I],D1,m[D1][i-1][j-1]+gcp[J])+sub;
	  t[M1][i][j]=trace;
	  
	  
	  m[D1][i][j]=dp_max (&trace,2, M1,m[M1][i][j-1]+gop[J]+gep[J],D1, m[D1][i][j-1]+gep[J]);
	  t[D1][i][j]=trace;
	  
	  
	  m[I1][i][j]=dp_max (&trace,2, M1,m[M1][i-1][j]+gop[I]+gep[I],I1, m[I1][i-1][j]+gep[I]);
	  t[I1][i][j]=trace;
	  
	}
	  
    }
  score=dp_max (&trace,3, M1,m[M1][li][lj],D1,m[D1][li][lj],I1, m[I1][li][lj]);
  
  LEN=0;i=li;j=lj;
  al=declare_char (2, li+lj);
  

  trace=t[trace][i][j];
  while (!(i==0 &&j==0))
    {
  
      ntrace=t[trace][i][j];
     
      
      if (i==0)
	{
	  al[0][LEN]=0;
	  al[1][LEN]=1;
	  j--;
	  LEN++;
	}
      else if ( j==0)
	{
	  al[0][LEN]=1;
	  al[1][LEN]=0;
	  i--;
	  LEN++;
	}
      else if ( trace==M1)
	{
	  al[0][LEN]=1;
	  al[1][LEN]=1;
	  i--; j--;
	  LEN++;
	}
      else if ( trace==D1)
	{
	  al[0][LEN]=0;
	  al[1][LEN]=1;
	  j--;
	  LEN++;
	}
      else if ( trace == I1)
	{
	  al[0][LEN]=1;
	  al[1][LEN]=0;
	  i--;
	  LEN++;
	}
      trace=ntrace;     
      
    }
  
  invert_list_char ( al[0], LEN);
  invert_list_char ( al[1], LEN);	
  if ( A->declared_len<=LEN)A=realloc_aln  ( A,2*LEN+1);	
  
  aln=A->seq_al;
  char_buf= vcalloc (LEN+1, sizeof (char));	
  for ( c=0; c< 2; c++)
    {
      for ( a=0; a< ns[c]; a++) 
	{		
	  int ch=0;
	  for ( b=0; b< LEN; b++)
	    {		   
	      if (al[c][b]==1)
		char_buf[b]=aln[l_s[c][a]][ch++];
	      else
		char_buf[b]='-';
	    }
	  char_buf[b]='\0';
	  sprintf (aln[l_s[c][a]],"%s", char_buf);
	}
    }
  
  
  A->len_aln=LEN;
  A->nseq=ns[0]+ns[1];
  vfree (char_buf);
  free_char (al, -1);
  free_int (pos0, -1);
  return score;
}
/*******************************************************************************/
/*                GLOCAL 2                                                     */
/*                                                                             */
/*	makes DP between the the ns[0] sequences and the ns[1]                 */
/*                                                                             */
/*	for MODE, see the function  get_dp_cost                                */
/*******************************************************************************/
int glocal2_pair_wise (Alignment *IN,int*ns, int **ls,Constraint_list *CL)
{
  int a, b, s=0;
  Alignment *A, *R,*L;
  char *seq, *buf;
  
  buf=vcalloc (1000, sizeof (char));
  seq=vcalloc (1000, sizeof (char));
  
  A=copy_aln (IN,NULL);
  L=copy_aln (IN,NULL);
  R=copy_aln (IN,NULL);
  
  gotoh_pair_wise_sw (A, ns, ls, CL);
  
  HERE ("1");
  for (a=0; a<2; a++)
    {
      for (b=0; b<ns[a]; b++)
	{
	  s=ls[a][b];
	  sprintf ( seq,"%s", IN->seq_al[s]);
	  
	  seq[A->order[s][2]]='\0';
	  sprintf (L->seq_al[s], "%s", seq);
	  sprintf (R->seq_al[s], "%s", seq+A->order[s][3]+1);
	}
    }
  HERE ("2");
  print_sub_aln (A, ns, ls);
  gotoh_pair_wise(L, ns, ls, CL);
  print_sub_aln (L, ns, ls);
  gotoh_pair_wise(R, ns, ls, CL);
  print_sub_aln (R, ns, ls);
  
  IN=realloc_aln (IN, A->len_aln+L->len_aln+R->len_aln+1);
  for (a=0; a<2; a++)
    {
      for (b=0; b<ns[a]; b++)
	{
	  s=ls[a][b];
	  sprintf (IN->seq_al[s], "%s%s%s",L->seq_al[s], A->seq_al[s], R->seq_al[s]);
	}
    }
  IN->len_aln=strlen (IN->seq_al[s]);
  
  print_sub_aln (IN, ns, ls);
  vfree (seq); vfree (buf);
  free_aln (A); free_aln (L);free_aln (R);
  return IN->score_aln;
}


int gotoh_pair_wise (Alignment *A,int*ns, int **l_s,Constraint_list *CL)
	{
/*******************************************************************************/
/*                NEEDLEMAN AND WUNSCH (GOTOH)                                 */
/*                                                                             */
/*	makes DP between the the ns[0] sequences and the ns[1]                 */
/*                                                                             */
/*	for MODE, see the function  get_dp_cost                                */
/*******************************************************************************/
	  

/*TREATMENT OF THE TERMINAL GAP PENALTIES*/
/*TG_MODE=0---> gop and gep*/
/*TG_MODE=1---> ---     gep*/
/*TG_MODE=2---> ---     ---*/


	    int TG_MODE;
	    int l_gop, l_gep;
	    int gop, gep;
	    int maximise;
/*VARIANLES FOR THE MULTIPLE SEQUENCE ALIGNMENT*/
	int a, b, i, j;

	int *cc;	
	int *dd,*ddg;
	int   e, eg;

	int lenal[2], len;
	int t, c=0,s, ch;
	int sub;
	int fop;
	int score=0;
        int **pos0;
	static char **al;
	char **aln;
	int ala, alb,LEN;
	char *buffer;
	char *char_buf;
/*trace back variables       */
	FILE       *long_trace=NULL;
	TRACE_TYPE *buf_trace=NULL;
	static TRACE_TYPE **trace;
	TRACE_TYPE k;
	TRACE_TYPE *tr;
	int long_trace_flag=0;
	int dim;
/********Prepare penalties*******/
	gop=CL->gop*SCORE_K;
	gep=CL->gep*SCORE_K;
	TG_MODE=CL->TG_MODE;
	maximise=CL->maximise;
	
	
/********************************/	
/*CLEAN UP AFTER USE*/
	if ( A==NULL)
	   {
	   free_int (trace,-1);
	   trace=NULL;
	   free_char (al,-1);
	   al=NULL;
	   return 0;
	   }		

/*DO MEMORY ALLOCATION FOR DP*/

	lenal[0]=strlen (A->seq_al[l_s[0][0]]);
	lenal[1]=strlen (A->seq_al[l_s[1][0]]);
	len= MAX(lenal[0],lenal[1])+1;
	buf_trace=vcalloc ( len, sizeof (TRACE_TYPE));	
	buffer=vcalloc ( 2*len, sizeof (char));	
        al=declare_char (2, 2*len);  
	
	char_buf= vcalloc (2*len, sizeof (char));	
	

	dd = vcalloc (len, sizeof (int));
	

	cc = vcalloc (len, sizeof (int));
	ddg=vcalloc (len, sizeof (int));
	

	
	if ( len>=MAX_LEN_FOR_DP)
	    {
	    long_trace_flag=1;
	    long_trace=vtmpfile();
	    }
	else
	    {
	   
	    dim=(trace==NULL)?0:read_size_int ( trace,sizeof (int*));	   
	    trace    =realloc_int ( trace,dim,dim,MAX(0,len-dim), MAX(0,len-dim));
	    }
	
/*END OF MEMORY ALLOCATION*/
	
	
		/*
		0(s)   +(dd)
		  \      |
		   \     |
		    \    |
		     \   |
		      \  |
		       \ |
		        \|
		-(e)----O
		*/ 
		       
	pos0=aln2pos_simple ( A,-1, ns, l_s);


	cc[0]=0;		
	tr=(long_trace_flag)?buf_trace:trace[0];
	tr[0]=(TRACE_TYPE)1;
	for ( j=1; j<=lenal[1]; j++)tr[j]=(TRACE_TYPE)-1;
	if (long_trace_flag)fwrite (buf_trace, sizeof ( TRACE_TYPE),lenal[1]+1, long_trace);
	
	
	t=(TG_MODE==0)?gop:0;
	

	for (cc[0]=0,j=1; j<=lenal[1]; j++)
	    {
	    
	    l_gop=(TG_MODE==0)?gop:0;
	    l_gep=(TG_MODE==2)?0:gep;

	    cc[j]=t=t+l_gep;
	    dd[j]=  t+  gop;
	    }

	t=(TG_MODE==0)?gop:0;	
	
	for (i=1; i<=lenal[0];i++)
			{			
			tr=(long_trace_flag)?buf_trace:trace[i];
			s=cc[0];

			l_gop=(TG_MODE==0)?gop:0;
			l_gep=(TG_MODE==2)?0:gep;
			
			
			
			cc[0]=c=t=t+l_gep;
			e=t+  gop;
			tr[0]=(TRACE_TYPE)1;

			

			for (eg=0,j=1; j<=lenal[1];j++)
				{				   
				 
				  sub=(CL->get_dp_cost) (A, pos0, ns[0], l_s[0], i-1, pos0, ns[1], l_s[1],j-1,CL);	
				      
				/*get the best Insertion*/
				l_gop=(i==lenal[0] || i==1 )?((TG_MODE==0)?gop:0):gop;
				l_gep=(i==lenal[0] || i==1)?((TG_MODE==2)?0:gep):gep;
			

				if ( a_better_than_b ( e,c+l_gop, maximise))eg++;
				else eg=1;	
				e=best_of_a_b (e, c+l_gop, maximise)+l_gep;
				
				/*Get the best deletion*/
				l_gop=(j==lenal[1] || j==1)?((TG_MODE==0)?gop:0):gop;
				l_gep=(j==lenal[1] || j==1)?((TG_MODE==2)?0:gep):gep;
				

				if ( a_better_than_b ( dd[j], cc[j]+l_gop, maximise))ddg[j]++;
				else ddg[j]=1;
				dd[j]=best_of_a_b( dd[j], cc[j]+l_gop,maximise)+l_gep;
				


				c=best_int(3,maximise,&fop, e, s+sub,dd[j]);
				/*Chose Substitution for tie breaking*/
				if ( fop==0 && (s+sub)==e)fop=1;
				else if ( fop==2 && (s+sub)==dd[j])fop=1;
				/*Chose Deletion for tie breaking*/
				else if ( fop==2 && e==dd[j])fop=1;

				fop-=1;
				s=cc[j];
				cc[j]=c;	

	
				if ( fop<0)
					{tr[j]=(TRACE_TYPE)fop*eg;
					}
				else if ( fop>0)
				        {tr[j]=(TRACE_TYPE)fop*ddg[j];
					}
				else if (fop==0)
					{tr[j]=(TRACE_TYPE)0;	
					}					
				fop= -2;
				}
			if (long_trace_flag)
			    {
			    fwrite ( buf_trace, sizeof (TRACE_TYPE), lenal[1]+1, long_trace);
			    }
			}
	
	score=c;
	
        i=lenal[0];
	j=lenal[1];
	ala=alb=0;
	

	while (i>=0 && j>=0 && ((i+j)!=0))
			{
			if ( i==0)
				k=-1;
			else if ( j==0)
				k=1;
			else if ( j==0 && i==0)
				k=1;	
			else
			        {
				if (long_trace_flag)
				   {
				   fseek ( long_trace, sizeof (TRACE_TYPE)*((lenal[1]+1)*(i)+j),SEEK_SET);
				   fread ( &k, sizeof (TRACE_TYPE), 1, long_trace);
				   }
				else
				   {
				   
				   k=trace[i][j];
				   }
				}
				
				
			if (k==0)
				{
				
				al[0][ala++]=1;
				al[1][alb++]=1;
				i--;
				j--;
				}		
			else if (k>0)
				{
				
				for ( a=0; a< k; a++)
					{
					al[0][ala++]=1;
					al[1][alb++]=0;
					i--;
					}
				}
			else if (k<0)
				{
				
				for ( a=0; a>k; a--)
					{
					al[0][ala++]=0;
					al[1][alb++]=1;
					j--;
					}
				}
			}
      
	LEN=ala;	
	c=LEN-1;  
	
	

	invert_list_char ( al[0], LEN);
	invert_list_char ( al[1], LEN);	
	if ( A->declared_len<=LEN)A=realloc_aln2  ( A,A->max_n_seq, 2*LEN);	
	aln=A->seq_al;

	for ( c=0; c< 2; c++)
	    {
	    for ( a=0; a< ns[c]; a++) 
		{		
		ch=0;
		for ( b=0; b< LEN; b++)
		    {		   
		    if (al[c][b]==1)
			char_buf[b]=aln[l_s[c][a]][ch++];
		    else
			char_buf[b]='-';
		   }
		char_buf[b]='\0';
		sprintf (aln[l_s[c][a]],"%s", char_buf);
	        }
	     }
	
	
	A->len_aln=LEN;
	A->nseq=ns[0]+ns[1];
	

	vfree ( cc);
	vfree (dd);		
	vfree (ddg);
	vfree (buffer);
	vfree (char_buf); 
	vfree (buf_trace);
	free_char ( al, -1);
	free_int (pos0, -1);
	if ( long_trace_flag)fclose (long_trace);	



	return score;
	}
     

int get_transition_cost (Alignment *A, int **posi, int ni, int *li, int i, int **posj, int nj, int *lj, int j,Constraint_list *CL);
int gotoh_pair_wise_lgp_sticky ( Alignment *A, int *ns, int **l_s, Constraint_list *CL)
{
  int i,j, li, lj, n, sub, trace,ntrace, a, b, c, score;
  int I, J;
  int M1, I1, D1, LEN;
  char **al, *char_buf, **aln;
  int **pos0;
  
  int gop[2], gcp[2], gep[2];
  static int ***gpl, ***t, ***m;
  static int max_li, max_lj;
  
  

  //gotoh_pair_wise ( A, ns, l_s,CL);
  //ungap_sub_aln (A, ns[0], l_s[0]);
  //ungap_sub_aln (A, ns[1], l_s[1]);
  
  I=0;J=1;

  
  li=strlen (A->seq_al[l_s[I][0]]);
  lj=strlen (A->seq_al[l_s[J][0]]);
  
  if ( !gpl)gpl=vcalloc ( 2, sizeof (int**));
  gpl[I]=aln2local_penalties (A,ns[I], l_s[I], CL,gpl[I]);
  gpl[J]=aln2local_penalties (A,ns[J], l_s[J], CL,gpl[J]);
  
  
  n=1;
  M1=n++;D1=n++;I1=n++;
  
  if ( li>max_li ||lj>max_lj )
    {
      free_arrayN((void**)t, 3);
      free_arrayN((void**)m, 3);

     
      max_li=li;
      max_lj=lj;
      t=declare_arrayN(3, sizeof (int),n, max_li+1, max_lj+1);
      m=declare_arrayN(3, sizeof (int),n, max_li+1, max_lj+1);
      
    }
  pos0=aln2pos_simple ( A,-1, ns, l_s);
 
  
  for (j=1; j<=lj; j++)
    {
      gep[J]=gpl[J][GEP][j-1];
      m[D1][0][j]=gep[J]*j;
      m[I1][0][j]=m[D1][0][j]-1;
      m[M1][0][j]=m[D1][0][j]-1;
    }
  
  //D1: gap in sequence I
  //I1: gap in sequence J
  
  
  for (i=1; i<=li; i++)
    {
      gep[I]=gpl[I][GEP][i-1];
      gop[I]=gpl[I][GOP][i-1];
      gcp[I]=gpl[I][GCP][i-1];
      
      m[I1][i][0]=i*gep[I];
      m[D1][i][0]= m[I1][i][0]-1;
      m[M1][i][0]= m[I1][i][0]-1;
                 
     
      
      gop[I]=(i==1 || i==li )?0:gop[I];
      gcp[I]=(i==1 || i==li )?0:gcp[I];
      
      
      for ( j=1; j<=lj; j++)
	{
	  int transition;
	  
	  gep[J]=gpl[J][GEP][j-1];
	  gop[J]=gpl[J][GOP][j-1];
	  gcp[J]=gpl[J][GCP][j-1];
	  
	  //gep[J]=gep[I]=(gep[J]+gep[I])/2;
	  //gop[J]=gop[I]=(gop[J]+gop[I])/2;
	  //gcp[J]=gcp[I]=(gcp[J]+gcp[I])/2;
	  

	  gop[J]=(j==1 || j==lj )?0:gop[J];
	  gcp[J]=(j==1 || j==lj )?0:gcp[J];
	  
	  
	  sub=(CL->get_dp_cost) (A, pos0, ns[0], l_s[0], i-1, pos0, ns[1], l_s[1],j-1,CL);
	  transition=get_transition_cost (A, pos0, ns[0], l_s[0], i-1, pos0, ns[1], l_s[1],j-1,CL);

	  m[M1][i][j]=dp_max (&trace,3,M1,m[M1][i-1][j-1]+transition,I1, m[I1][i-1][j-1]+gcp[I],D1,m[D1][i-1][j-1]+gcp[J])+sub;
	  t[M1][i][j]=trace;
	  
	  
	  m[D1][i][j]=dp_max (&trace,2, M1,m[M1][i][j-1]+gop[J]+gep[J],D1, m[D1][i][j-1]+gep[J]);
	  t[D1][i][j]=trace;
	  
	  
	  m[I1][i][j]=dp_max (&trace,2, M1,m[M1][i-1][j]+gop[I]+gep[I],I1, m[I1][i-1][j]+gep[I]);
	  t[I1][i][j]=trace;
	  
	}
	  
    }
  score=dp_max (&trace,3, M1,m[M1][li][lj],D1,m[D1][li][lj],I1, m[I1][li][lj]);
  
  LEN=0;i=li;j=lj;
  al=declare_char (2, li+lj);
  

  trace=t[trace][i][j];
  while (!(i==0 &&j==0))
    {
  
      ntrace=t[trace][i][j];
     
      
      if (i==0)
	{
	  al[0][LEN]=0;
	  al[1][LEN]=1;
	  j--;
	  LEN++;
	}
      else if ( j==0)
	{
	  al[0][LEN]=1;
	  al[1][LEN]=0;
	  i--;
	  LEN++;
	}
      else if ( trace==M1)
	{
	  al[0][LEN]=1;
	  al[1][LEN]=1;
	  i--; j--;
	  LEN++;
	}
      else if ( trace==D1)
	{
	  al[0][LEN]=0;
	  al[1][LEN]=1;
	  j--;
	  LEN++;
	}
      else if ( trace == I1)
	{
	  al[0][LEN]=1;
	  al[1][LEN]=0;
	  i--;
	  LEN++;
	}
      trace=ntrace;     
      
    }
  
  invert_list_char ( al[0], LEN);
  invert_list_char ( al[1], LEN);	
  if ( A->declared_len<=LEN)A=realloc_aln  ( A,2*LEN+1);	
  
  aln=A->seq_al;
  char_buf= vcalloc (LEN+1, sizeof (char));	
  for ( c=0; c< 2; c++)
    {
      for ( a=0; a< ns[c]; a++) 
	{		
	  int ch=0;
	  for ( b=0; b< LEN; b++)
	    {		   
	      if (al[c][b]==1)
		char_buf[b]=aln[l_s[c][a]][ch++];
	      else
		char_buf[b]='-';
	    }
	  char_buf[b]='\0';
	  sprintf (aln[l_s[c][a]],"%s", char_buf);
	}
    }
  
  
  A->len_aln=LEN;
  A->nseq=ns[0]+ns[1];
  vfree (char_buf);
  free_char (al, -1);
  free_int (pos0, -1);
  return score;
}
int get_transition_cost (Alignment *A, int **posi, int ni, int *li, int i, int **posj, int nj, int *lj, int j,Constraint_list *CL)
{
  /*counts the number of identical transitions between position i-1, i and j-1..j*/
  float t=0;
  int a,s;
  Sequence *S;
  
  if (i==0 || j==0)return 0;
  
  for (a=0; a<ni; a++)
    {
      s=li[a];
      if (posi[s][i]<0 || posi[s][i-1]<0)continue;
      if (S->seq[li[a]][i-1]==S->seq[li[a]][i-1])t++;
    }
  
  for (a=0; a<nj; a++)
    {
      s=lj[a];
      if (posj[s][j]<0 || posj[s][j-1]<0)continue;
      if (S->seq[li[a]][j-1]==S->seq[li[a]][j-1])t++;
    }

  t=(t*10)/(float)(ni+nj);
  return t;
}
/*******************************************************************************/
/*                idscore_pairseq: measure the % id without delivering thze aln*/                                                   
/*                                                                             */
/*	makes DP between the the ns[0] sequences and the ns[1]                 */
/*                                                                             */
/*	for MODE, see the function  get_dp_cost                                */
/*******************************************************************************/
int cl2pair_list ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in, int mode, int ndiag);
int cl2pair_list_ref ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in);
int cl2pair_list_ecf ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in);
int cl2pair_list_diag ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in, int add);
int cl2list_borders   (Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in);
int cl2diag_cap (Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in); //add one element at the end of each segment so that they can be joined
int** cl2sorted_diagonals   ( Alignment *A, int *ns, int **ls, Constraint_list *CL);
int** cl2sorted_diagonals_mat  ( Alignment *A, int *ns, int **ls, Constraint_list *CL);
int** cl2sorted_diagonals_cs   ( Alignment *A, int *ns, int **ls, Constraint_list *CL);
int list2nodup_list (Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in);
int fill_matrix ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in);
int cl2pair_list ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in, int mode, int ndiag)
{
  int v;
  if (!A)
    {
      free_int (list_in[0], -1); list_in[0]=NULL;
      n_in[0]=0;
    }
  
  cl2list_borders(A, ns, ls, CL, list_in, n_in);
  
  if ( mode==0)
    v=cl2pair_list_ref (A, ns, ls, CL, list_in, n_in);
  else if (mode==1)
    v=cl2pair_list_ecl (A, ns, ls, CL, list_in, n_in);
  else if (mode==2)
    v=cl2pair_list_diag (A, ns, ls, CL, list_in, n_in,ndiag); //add diagonals
  
  cl2diag_cap (A, ns, ls, CL, list_in, n_in);
  //fill_matrix (A, ns, ls, CL, list_in, n_in);//Fill matrix with 0s
  sort_list_int (list_in[0],7, 1, 0, n_in[0]-1);
  list2nodup_list (A, ns, ls, CL, list_in, n_in);
  return v;
  
}
int fill_matrix( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list, int *n)
{
  int a, b, l1, l2, n2=0;
  int score;
  int **pos;
  int max_n;
  if (!A) return 0;
  pos=aln2pos_simple ( A,-1, ns, ls);
  l1=strlen (A->seq_al[ls[0][0]]);
  l2=strlen (A->seq_al[ls[1][0]]);
  
  max_n=read_array_size (list[0], sizeof (int));
  
  for (a=0; a<=l1; a++)
    for (b=0; b<=l2; b++)
      {
	score=0;
	score=(a==0 || b==0)?0:slow_get_dp_cost (A, pos, ns[0], ls[0],a-1, pos, ns[1], ls[1], b-1, CL);
	if ( score>0 && a!=0 && b!=0 && a!=l1 && b!=l2)
	  {
	    if (n[0]==max_n){max_n+=1000;list[0]=vrealloc (list[0], max_n*sizeof (int*));}
	    if (!list[0][n[0]])list[0][n[0]]=vcalloc (7, sizeof (int));
	    list[0][n[0]][0]=a;
	    list[0][n[0]][1]=b;
	    list[0][n[0]][3]=(l1-a)+b;
	    list[0][n[0]][2]=score;
	    if ( a!=0 && b!=0 && a!=l1 && b!=l2)
	      {
		n2++;
	      }
	    n[0]++;
	  }

      }

  return n[0];
  }
int list2nodup_list ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in)
{
  int **list;
  int n, a, b, c;
  
  list=list_in[0];
  n=n_in[0];
  
  if ( !A)return 0;
  
  for (b=a=1; a<n; a++)
    {
      if (list[a][0]==list[b-1][0] && list[a][1]==list[b-1][1])
	{
	  //HERE ("Duplicate");
	   list[b-1][2]=MAX(list[b-1][2],list[a][2]);
	}
      else
	{
	  for (c=0; c<4; c++)list[b][c]=list[a][c];
	  b++;
	}
	
	}
  n_in[0]=b;
  list_in[0]=list;
  return b;
}
int** cl2sorted_diagonals   ( Alignment *A, int *ns, int **ls, Constraint_list *CL)
{
  if ( CL && CL->L)return cl2sorted_diagonals_cs   ( A, ns, ls, CL);
  else return cl2sorted_diagonals_mat   ( A, ns, ls, CL);
}

static int kword;
static char **warray;
int cmp_word ( const int**a, const int**b);
int ** seq2index_list ( Sequence *S, int k);
int** cl2sorted_diagonals_mat  ( Alignment *A, int *ns, int **ls, Constraint_list *CL)
{
  
  int a,b,c,d, comp, k, l1, l2, ndiag;
  int  **diag;
  
  
  static char **alp, alps=5;
  char *buf1, *buf2;
  
  if (!A)return NULL;
  if ( !alp)
    alp=make_group_aa_upgma ("blosum62mt",alps);
  
  
  k=2;
  l1=strlen (A->seq_al[ls[0][0]]);
  buf1=vcalloc ( l1+1, sizeof (char));
  l2=strlen (A->seq_al[ls[1][0]]);
  buf2=vcalloc ( l2+1, sizeof (char));
  
  ndiag=l1+l2;
  diag=declare_int (ndiag+3,2);
  for (a=0; a<=ndiag; a++) diag[a][0]=a;
  vfree (diag[ndiag+1]);
  diag[ndiag+1]=NULL;
 
  for ( a=0; a<ns[0]; a++)
    {
      sprintf (buf1, A->seq_al[ls[0][a]]);
      lower_string (buf1);
      string_convert (buf1, alps, alp);
      for (b=0; b<ns[1]; b++)
	{
	  sprintf (buf2, A->seq_al[ls[1][b]]);
	  lower_string (buf2);
	  string_convert (buf2, alps, alp);
	  for (c=0; c<l1-k; c++)
	    {
	      if ( strnrchr (buf1+c, '-', k))continue;
	      for (d=0; d<l2-k; d++)
	      {
		if ( strnrchr (buf2+d, '-', k))continue;
		comp=strncmp (buf1+c,buf2+d, k);
		diag[l1-c+d][1]+=(comp==0)?1:0;
	      }
	    }
	}
    }
  HERE ("dONE");
  /*
   max_len=MAX(l1, l2);
   for (a=1; a<ndiag; a++)
    {
      int l, d, p1_0, p1_1, p2_0, p2_1;
      d=diag[a][0];

      if (d<=l1)l=MIN(d,l2);
      else l=MIN(((l1+l2)-d),l1); 

      diag[a][1]=(diag[a][1]*1000)/max_len;
      diag[a][1]=(float)((float)diag[a][1]*(((float)max_len)/(float)l));
    }
  */
  sort_int_inv (diag, 2, 1,0,ndiag);
  
  return diag;
}







int ** seq2index_list ( Sequence *S, int k)
{
  int **list,**mlist=NULL;
  int a, b,c,ml, n, l, e, s, max=0, nm=0;
  char *cw;
  

  for (ml=0,a=0; a<S->nseq; a++)ml+=strlen (S->seq[a]);
  list=declare_int (ml+1, 2);
  
  for (n=0,a=0; a<S->nseq; a++)
    {
      l=strlen (S->seq[a])-k;
      for ( b=0; b<l; b++, n++)
	{
	  list[n][0]=a;
	  list[n][1]=b;
	}
    }
  list[n][0]=-1;
  
  warray=S->seq;
  kword=k;
  qsort (list, n, sizeof (long**), (int(*)(const void*,const void*))(cmp_word));
    
  cw=NULL;
  e=s=0;
  nm=0;
  for (a=0; a<=n; a++)
    {
      int s1, s2, r1, r2;
      if (!cw ||a==n|| strncmp (warray[list[a][0]]+list[a][1],cw, k)!=0)
	{
	  if (a<n)cw=warray[list[a][0]]+list[a][1];
	  for (b=s; b<a-1; b++)
	    for (c=b+1; c<a; c++)
	      {

		if (list[b][0]<list[c][0])
		  {s1=list[b][0];
		    r1=list[b][1];
		    s2=list[c][0];
		    r2=list[c][1];
		  }
		else
		  {
		    s2=list[b][0];
		    r2=list[b][1];
		    s1=list[c][0];
		    r1=list[c][1];
		  }
		
		if (s1==s2)continue;
		else
		  {
		    if (nm>=max){max+=1000; mlist=vrealloc (mlist, max*sizeof (int*));}
		    mlist[nm]=vcalloc (4, sizeof (int));
		    mlist[nm][0]=s1;
		    mlist[nm][1]=s2;
		    mlist[nm][2]=r1;
		    mlist[nm][3]=r2;
		    nm++;
		  }
		s=a;
	      }
	}
    }
  
  if (nm>=max){max+=1000;mlist=vrealloc (mlist, max*sizeof (int));}
  sort_list_int ( mlist,4,1, 0, nm-1);
  return mlist;
}
int cmp_word ( const int**a, const int**b)
{
  int c;

  c=strncmp (warray[a[0][0]]+a[0][1], warray[b[0][0]]+b[0][1], kword);

  
  if (c) return c;
  else 
    {
      for (c=0; c<2; c++)
	{
	  if ( a[0][c]>b[0][c])return   1;
	  else if (a[0][c]<b[0][c])return   -1;
	}
    }
  return 0;
}
 
int** cl2sorted_diagonals_cs   ( Alignment *A, int *ns, int **ls, Constraint_list *CL)
{
  int p1, p2, s1, s2, r1, r2;
  int a, b, l1, l2;
  int **pos;
  int **diag;
  int ndiag;
  int diag_i;

  
  pos=aln2pos_simple ( A,-1, ns, ls);
  
  l1=strlen (A->seq_al[ls[0][0]]);
  l2=strlen (A->seq_al[ls[1][0]]);
  
  

  CL=index_res_constraint_list (CL, CL->weight_field);
  ndiag=l1+l2;
  diag=declare_int (ndiag+3, 2);
  
  for (a=1; a<=ndiag; a++)diag[a][0]=a;
    
  for (p1=1; p1<=l1; p1++)
    {
      for (p2=1; p2<=l2; p2++)
	{
	  for (a=0; a<ns[0]; a++)
	    {
	      s1=ls[0][a];
	      r1=pos[s1][p1-1];
	      if (r1<=0)continue;
	      for (b=0; b<ns[1]; b++)
		{
		  s2=ls[1][b];
		  r2=pos[s2][p2-1];
		  if (r2<=0)continue;
		  diag_i=(l1-p1)+p2;
		  
		  diag[diag_i][1]+=residue_pair_extended_list_raw (CL,s1, r1-1,s2, r2-1);
		}
	    }
	}
    }
  
  sort_int_inv (diag, 2, 1,0,ndiag);

  vfree (diag[ndiag+1]);
  diag[ndiag+1]=NULL;
  free_int (pos, -1);
  return diag;
} 
int** cl2sorted_diagonals_cs_old2   ( Alignment *A, int *ns, int **ls, Constraint_list *CL)
{
  int p1, p2, si, s, r, t_s, t_r;
  int a, l1, l2;
  int *sl2, **pos,**inv_pos;
  int **diag;
  int ndiag;
  int diag_i;
  int max_len;
  
  pos=aln2pos_simple ( A,-1, ns, ls);
  inv_pos=vcalloc ((CL->S)->nseq, sizeof (int*));
  for (a=0; a<ns[1]; a++)inv_pos[ls[1][a]] =seq2inv_pos(A->seq_al[ls[1][a]]);
  
  l1=strlen (A->seq_al[ls[0][0]]);
  l2=strlen (A->seq_al[ls[1][0]]);
  sl2=vcalloc ((CL->S)->nseq, sizeof (int));
  
  for (a=0;a<ns[1]; a++)sl2[ls[1][a]]=1;
  CL=index_res_constraint_list (CL, CL->weight_field);
  ndiag=l1+l2;
  diag=declare_int (ndiag+3, 2);
  
  for (a=1; a<=ndiag; a++)diag[a][0]=a;
  for (p1=0; p1<=l1; p1++)
    {
      for (si=0;p1>0 && si<ns[0]; si++)
	{
	  s=ls [0][si];
	  r=pos[s][p1-1];
	  
	  for (a=1; r>0 && a<CL->residue_index[s][r][0];a+=3)
	    {
	      
	      t_s=CL->residue_index[s][r][a];
	      t_r=CL->residue_index[s][r][a+1];
	      
	      if (sl2[t_s])
		{
		  p2=inv_pos[t_s][t_r];
		  diag_i=(l1-p1)+p2;
		  diag[diag_i][1]+=(CL->get_dp_cost) (A, pos, ns[0], ls[0], p1-1, pos, ns[1], ls[1],p2-1,CL);	
		}
	    }
	}
    }
  max_len=MAX(l1, l2);
  for (a=1; a<ndiag; a++)
    {
      int l, d;
      d=diag[a][0];

      if (d<=l1)l=MIN(d,l2);
      else l=MIN(((l1+l2)-d),l1); 

      diag[a][1]/=max_len;
      diag[a][1]=(float)((float)diag[a][1]*(((float)max_len)/(float)l));
    }
  sort_int_inv (diag, 2, 1,0,ndiag);

  vfree (diag[ndiag+1]);
  diag[ndiag+1]=NULL;
  free_int (pos, -1);
  free_int (inv_pos, -1);
  vfree (sl2);
  
  return diag;
}
int cl2list_borders   (Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in)
{
  int maxlen, a,n, p1, p2, l1, l2;
  int **list;
  int **pos;
  if (!A)return 0;
  
//   printf("ICH BIN HIER\n");
  if ( list_in[0] && list_in[0][0]==0)
	  return read_array_size (list, sizeof (int*));
  
//   printf("UND HIER IRGENDWIE AUCH\n");
  
  list=list_in[0];
  n=n_in[0];
  if (!list)maxlen=0;
  else maxlen=read_array_size (list, sizeof (int*));
  

  l1=strlen (A->seq_al[ls[0][0]]);
  l2=strlen (A->seq_al[ls[1][0]]);
//    pos=aln2pos_simple ( A,-1, ns, ls);

  for (p1=0; p1<=l1; p1++)
    {
      if (p1==0 || p1==l1)
	{
	  for (p2=0; p2<=l2; p2++)
	    {
	      if (n==maxlen){maxlen+=1000;list=vrealloc (list,maxlen*sizeof (int*));}
	      if (!list[n])list[n]=vcalloc (7, sizeof (int));
	      list[n][0]=p1;
	      list[n][1]=p2;
	      list[n][3]=(l1-(p1))+(p2);
	      //list[n][2]=(p1==0||p2==0)?0:(CL->get_dp_cost) (A, pos, ns[0], ls[0], p1-1, pos, ns[1], ls[1],p2-1,CL);;
	      list[n][2]=(CL->gep)*SCORE_K*p2;
	      n++;
	    }
	}
      else
	{
	  for (a=0; a<2; a++)
	    {
	      p2=(a==0)?0:l2;
	      if (n==maxlen){maxlen+=1000;list=vrealloc (list,maxlen*sizeof (int*));}
	      if (!list[n])list[n]=vcalloc (7, sizeof (int));
	      list[n][0]=p1;
	      list[n][1]=p2;
	      list[n][3]=(l1-(p1))+(p2);
	      //list[n][2]=(p1==0||p2==0)?0:(CL->get_dp_cost) (A, pos, ns[0], ls[0], p1-1, pos, ns[1], ls[1],p2-1,CL);;
	      list[n][2]=(CL->gep)*SCORE_K*p1;
	      n++;
	    }
	}
    }
//   free_int (pos, -1);
  list_in[0]=list;
  n_in[0]=n;
  return read_array_size (list, sizeof (int*));
}

int cl2diag_cap (Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in)
{
  int **list;
  int n, in, a, b, al1, al2;
  int max_n;
  int cap=0;
  
  if (!A) return 0;
  
  al1=strlen (A->seq_al[ls[0][0]]);
  al2=strlen (A->seq_al[ls[1][0]]);
  
  list=list_in[0];
  n=n_in[0];
  max_n=read_array_size (list, sizeof (int*));
  
  
  
  
  for (a=0; a< n; a++)
    {
      b=list[a][3];
      list[a][3]=list[a][0];
      list[a][0]=b;
      
    }
  sort_list_int (list, 4, 1, 0, n-1);
  for (a=0; a< n; a++)
    {
      b=list[a][3];
      list[a][3]=list[a][0];
      list[a][0]=b;
    }
  

  in=n;
  
  for (a=0; a<in; a++)
    {
      int i, j, pi, pj, ni, nj;
      if (list[a][2]==0)continue;
      i=list[a][0];
      j=list[a][1];
      
      if (a==0){pi=-10;pj=-10;}
      else {pi=list[a-1][0];pj=list[a-1][1];}
      
      if (a==in-1){ni=-10; nj=-10;}
      else {ni=list[a+1][0]; nj=list[a+1][1];}
      
      
      if (i==0 || j==0);
      else if ( i==pi || j==pj);
      else if ( i-pi!=1 || j-pj!=1)
	{
	  
	   if (n==max_n){max_n+=1000;list=vrealloc (list, max_n*sizeof (int*));}
	   if (!list[n])list[n]=vcalloc (7, sizeof (int));
	   
	   list[n][0]=i-1;
	   list[n][1]=j-1;
	   list[n][3]=list[a][3];
	   list[n][2]=cap;
	   n++;
	}
    
    
      if (i==al1 || j==al2);
      else if ( i==ni || j==nj);
      else if ( ni-i!=1 || nj-j!=1)
	{
	  
	  if (n==max_n){max_n+=1000;list=vrealloc (list, max_n*sizeof (int*));}
	  if (!list[n])list[n]=vcalloc (7, sizeof (int));
	  
	  list[n][0]=i+1;
	  list[n][1]=j+1;
	  list[n][3]=list[a][3];
	  list[n][2]=cap;
	  n++;
	}
    
    }
  list_in[0]=list;
  n_in[0]=n;
  return n;
}
	  
int cl2pair_list_diag_mat ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in, int add );
int cl2pair_list_diag_cl ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in, int add );	  
int cl2pair_list_diag ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in, int add )
{
  if (CL->L)return cl2pair_list_diag_cl (A, ns, ls, CL, list_in, n_in, add);
  else return cl2pair_list_diag_mat (A, ns, ls, CL, list_in, n_in, add);
}
int cl2pair_list_diag_mat ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in, int add )
{
  int p1, p2, n,d;
  int a, l1, l2;
  int set=0;
  static int **pos;
  int max_n;
  static int **diag;
  int **list;

  if (A==NULL)
    {
      free_int (pos, -1);pos=NULL;
      free_int (diag, -1);diag=NULL;
      //free_int (list_in[0], -1); list_in[0]=NULL;
      return 0;
    }
  
  if ( !pos)
    {
      pos=aln2pos_simple ( A,-1, ns, ls);
      diag=cl2sorted_diagonals (A,ns,ls,CL);
    }
  
  list=list_in[0];
  n=n_in[0];
  max_n=read_array_size (list, sizeof (int**));
  
  
  
  l1=strlen (A->seq_al[ls[0][0]]);
  l2=strlen (A->seq_al[ls[1][0]]);
  
  d=0;
  if ( add)
    {
      while ( diag[d] && diag[d][1]==-1)d++;
      add+=d;
    }
  else
    {
      d=0;
      while (diag[add++]);
    }
  
  HERE ("Add %d diagonals, starts %d N=%d", add, d, n);
  
  for (d=0; d<add && diag[d]; d++)
    {
      int p1_0, p2_0;
      
      set=1;

      HERE ("\t S=%d", diag[d][1]);
      
      p1_0=MAX(0,l1-diag[d][0]);
      p2_0=MAX(0,diag[d][0]-l1);
      diag[d][1]=-1;
      
      for (p1=p1_0, p2=p2_0; p1<=l1 && p2<=l2; p1++,p2++)
	{
	  if (!BORDER(p1,l1, p2,l2) )
	    {

	      if (n==max_n){max_n+=1000;list=vrealloc (list, max_n*sizeof (int*));}
	      if (!list[n])list[n]=vcalloc (7, sizeof (int));
	      
	      list[n][0]=p1;
	      list[n][1]=p2;
	      list[n][3]=(l1-(p1))+(p2);
	      list[n][2]=(CL->get_dp_cost) (A, pos, ns[0], ls[0], p1-1, pos, ns[1], ls[1],p2-1,CL);
	      n++;
	    }
	}
    }
  HERE ("Addition Finished n=%d", n);
  if (!set) return 0;
  sort_list_int (list,4, 1, 0, n-1);
  
  list_in[0]=list;
  n_in[0]=n;
  HERE ("\nN=%d r=%.3f [l1=%d l2=%d]", n, (float)n/(float)(l1*l2), l1, l2);
  return max_n;
}

int cl2pair_list_diag_cl ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in, int add )
{
  int p1, p2,n,d;
  int l1, l2;
  int score, set=0;
  static int **pos;
  int max_n;
  static int **diag;
  int **list;

  if (A==NULL)
    {
      free_int (pos, -1);pos=NULL;
      free_int (diag, -1);diag=NULL;
      //free_int (list_in[0], -1); list_in[0]=NULL;
      return 0;
    }
  
  if ( !pos)
    {
      pos=aln2pos_simple ( A,-1, ns, ls);
      diag=cl2sorted_diagonals (A,ns,ls,CL);
    }
  
  list=list_in[0];
  n=n_in[0];
  max_n=read_array_size (list, sizeof (int**));
  
  
  
  l1=strlen (A->seq_al[ls[0][0]]);
  l2=strlen (A->seq_al[ls[1][0]]);
  if ( add==0)add=l1+l2;
  d=0;
  while ( diag[d] && diag[d][1]==-1)d++;
  HERE ("Add %d diagonals, starts %d N=%d", add, d, n);
  add+=d;
  for (d; d<add && diag[d]; d++)
    {
      int p1_0, p2_0;
      if (diag[d][1]==0)continue;
      set=1;

      HERE ("\t S=%d", diag[d][1]);
      
      p1_0=MAX(0,l1-diag[d][0]);
      p2_0=MAX(0,diag[d][0]-l1);
      diag[d][1]=-1;
      
      for (p1=p1_0, p2=p2_0; p1<=l1 && p2<=l2; p1++,p2++)
	{
	  if (!BORDER(p1,l1, p2,l2) && (score=(CL->get_dp_cost) (A, pos, ns[0], ls[0], p1-1, pos, ns[1], ls[1],p2-1,CL))!=0)
	    {

	      if (n==max_n){max_n+=1000;list=vrealloc (list, max_n*sizeof (int*));}
	      if (!list[n])list[n]=vcalloc (7, sizeof (int));
	      
	      list[n][0]=p1;
	      list[n][1]=p2;
	      list[n][3]=(l1-(p1))+(p2);
	      list[n][2]=score;
	      n++;
	    }
	}
    }
  HERE ("Addition Finished n=%d", n);
  if (!set) return 0;
  sort_list_int (list,4, 1, 0, n-1);
  
  list_in[0]=list;
  n_in[0]=n;
  HERE ("\nN=%d r=%.3f [l1=%d l2=%d]", n, (float)n/(float)(l1*l2), l1, l2);
  return max_n;
}
  
int cl2pair_list_ecl_norm ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in);
int cl2pair_list_ecl_raw ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in);
int cl2pair_list_ecl_pc ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in);
int cl2pair_list_ecl_noext_raw ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in);
int cl2pair_list_ecl_rna2 ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in);

int cl2pair_list_ecl_rawquad ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in);
int cl2pair_list_ecl ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in)
{
  int mode=5;
  
  if ( mode==1)return cl2pair_list_ecl_norm         (A, ns, ls, CL, list_in, n_in);
  else if ( mode==2)return cl2pair_list_ecl_raw     (A, ns, ls, CL, list_in, n_in);
  else if ( mode==3)return cl2pair_list_ecl_rawquad (A, ns, ls, CL, list_in, n_in);
  else if ( mode==4)return cl2pair_list_ecl_noext_raw (A, ns, ls, CL, list_in, n_in);
  else if ( mode==5)return cl2pair_list_ecl_pc     (A, ns, ls, CL, list_in, n_in);
}
int cl2pair_list_ecl_noext_raw ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in)
{
  int p1, p2, si, s, r, t_s2, t_r2, t_w2, n,n2;
  int a, b, l1, l2;
  int score;
  int **pos;
  int **list;
  int max_n;


  
  int nused;
  int *used_list, **used;
  int *sl2, **inv_pos;

  int filter=10;
  

  
  if ( !A) return 0;
  list=list_in[0];
  n=n_in[0];
  max_n=read_array_size (list, sizeof (int*));
  
 
  n2=0;
  pos=aln2pos_simple ( A,-1, ns, ls);
  inv_pos=vcalloc ((CL->S)->nseq, sizeof (int*));
  for (a=0; a<ns[1]; a++)inv_pos[ls[1][a]] =seq2inv_pos(A->seq_al[ls[1][a]]);
  
  l1=strlen (A->seq_al[ls[0][0]]);
  l2=strlen (A->seq_al[ls[1][0]]);
  sl2=vcalloc ((CL->S)->nseq, sizeof (int));
  
  for (a=0;a<ns[1]; a++)sl2[ls[1][a]]=1;
 
  CL=index_res_constraint_list (CL, CL->weight_field);
  
  used=declare_int (l2+1,2);
  used_list=vcalloc (l2+1, sizeof (int));
  nused=0;
    
  for (p1=0; p1<=l1; p1++)
    {
      for (nused=0,si=0;p1>0 && si<ns[0]; si++)
	{
	  s=ls [0][si];r=pos[s][p1-1];
	  for (a=1; r>0 && a<CL->residue_index[s][r][0];a+=3)
	    {
	      t_s2=CL->residue_index[s][r][a];t_r2=CL->residue_index[s][r][a+1];t_w2=CL->residue_index[s][r][a+2];
	      if (sl2[t_s2])
		{
		  p2=inv_pos[t_s2][t_r2];
		  score=t_w2;
		  if (!used[p2][1] && score>0)
		    {
		      used_list[nused++]=p2;
		    }
		  used[p2][0]+=score;
		  used[p2][1]++;
		}
	    }
	}
      for (a=0; a<nused; a++)
	{
	  
	  p2=used_list[a];
	 
	  score=used[p2][0]*SCORE_K;
	  used[p2][0]=used[p2][1]=0;
	  	  
	  if (n==max_n){max_n+=1000;list=vrealloc (list, max_n*sizeof (int*));}
	  if (score>filter && p1!=0 && p2!=0 && p1!=l1 && p2!=l2)
	    {
	      if (!list[n])list[n]=vcalloc (7, sizeof (int));
	      
	      list[n][0]=p1;
	      list[n][1]=p2;
	      list[n][3]=(l1-(p1))+(p2);
	      list[n][2]=score;
	      n++;
	    }
	}
    }
  

  vfree (used);
  vfree (used_list);
  free_int (inv_pos, -1);
  free_int (pos, -1);
  vfree (sl2);

  n_in[0]=n;
  list_in[0]=list;
  return 1;
}

int cl2pair_list_ecl_raw ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in)
{
  int p1, p2, si, s, r, t_s, t_r,t_w, t_s2, t_r2, t_w2, n,tot;
  int a, b, l1, l2;
  int **pos,**list;
  int max_n;
  
  
  int set, raw_max,nscore, score, nused;
  int *used_list, **used;
  int *sl2, **inv_pos;

  int filter1=0, filter2=0, max=0;
  int **nr;
  long tot_score=0, avg;
  int new_n=0;
  
  if ( !A) return 0;
  list=list_in[0];
  n=n_in[0];
  max_n=read_array_size (list, sizeof (int*));
  

  pos=aln2pos_simple ( A,-1, ns, ls);
  inv_pos=vcalloc ((CL->S)->nseq, sizeof (int*));
  for (a=0; a<ns[1]; a++)inv_pos[ls[1][a]] =seq2inv_pos(A->seq_al[ls[1][a]]);
  
  l1=strlen (A->seq_al[ls[0][0]]);
  l2=strlen (A->seq_al[ls[1][0]]);
  sl2=vcalloc ((CL->S)->nseq, sizeof (int));
  nr=declare_int (2, MAX(l1,l2)+1);
  
  for (a=0; a<l1;a++)
    for (b=0; b<ns[0]; b++)
      if (!is_gap(A->seq_al[ls[0][b]][a]))nr[0][a+1]++;
  for (a=0; a<l2;a++)
    for (b=0; b<ns[1]; b++)
      if (!is_gap(A->seq_al[ls[1][b]][a]))nr[1][a+1]++;
  
  for (a=0;a<ns[1]; a++)sl2[ls[1][a]]=1;
  
  CL=index_res_constraint_list (CL, CL->weight_field);
  
  used=declare_int (l2+1,2);
  used_list=vcalloc (l2+1, sizeof (int));
  nused=0;
    
  for (raw_max=0,p1=0; p1<=l1; p1++)
    {
      for (tot=0,nused=0,si=0;p1>0 && si<ns[0]; si++)
	{
	  s=ls [0][si];r=pos[s][p1-1];
	  for (a=1; r>0 && a<CL->residue_index[s][r][0];a+=3)
	    {
	      t_s=CL->residue_index[s][r][a];t_r=CL->residue_index[s][r][a+1];t_w=CL->residue_index[s][r][a+2];
	      for (b=0; b<CL->residue_index[t_s][t_r][0];)
		{
		  if (b==0){t_s2=t_s;t_r2=t_r;t_w2=t_w;b++;}
		  else 
		    { 
		      t_s2=CL->residue_index[t_s][t_r][b];t_r2=CL->residue_index[t_s][t_r][b+1];t_w2=CL->residue_index[t_s][t_r][b+2];b+=3;
		    }
		  
		  if (sl2[t_s2])
		    {
		      p2=inv_pos[t_s2][t_r2];
		      score=MIN(t_w,t_w2);
		      if (score<filter1)score=0;
		      //if ( score<3) continue;
		      if (!used[p2][1] && score>0)
			{
			  used_list[nused++]=p2;
			}
		      tot+=score;
		      used[p2][0]+=score;
		      used[p2][1]++;
		    }
		}
	    }
	}
      //set the threshold to 1/2 of the best normalised score
      
      for (filter2=0,set=0,a=0; a<nused; a++)
	{
	  
	  p2=used_list[a];
	  
	  score=used[p2][0];
	  nscore=(score*100)/tot;
	  if (set==0){filter2=nscore;set=1;}
	  filter2=MAX(nscore,filter2);
	  if ( score<0)HERE ("*********** %d", score);
	}
      filter2/=2;
      max+=nused;
      for (a=0; a<nused; a++)
	{
	  
	  p2=used_list[a];
	 
	  //score=used[p2][0];
	  nscore=(used[p2][0]*100)/tot; //Normalized score used for filtering
	  score =used[p2][0];
	  raw_max=MAX(score, raw_max);
	  
	  used[p2][0]=used[p2][1]=0;
	  
	  if (n==max_n){max_n+=1000;list=vrealloc (list, max_n*sizeof (int*));}
	  if (nscore>=filter2 && p1!=0 && p2!=0 && p1!=l1 && p2!=l2)
	    {
	      if (!list[n])list[n]=vcalloc (7, sizeof (int));
	      
	      list[n][0]=p1;
	      list[n][1]=p2;
	      list[n][3]=(l1-(p1))+(p2);
	      list[n][2]=score;
	      n++;
	      tot_score+=score;
	      new_n++;
	    }
	}
    }
  avg=tot_score/new_n;
  
  //CL->gop=-1*avg*3;CL->gep=0;
  HERE ("FILTER: %d->%d [THR=%d]", max, n-n_in[0], filter2);
 
  vfree (used);
  vfree (used_list);
  free_int (inv_pos, -1);
  free_int (pos, -1);
  vfree (sl2);

  n_in[0]=n;
  list_in[0]=list;
  return 1;
}


/**
 * Calculates scores for diagonal segments.
 * 
 * \param Alignment The sequences.
 * \param ns Number of sequences in each group
 * \param ls sequences in in groups (ls[0][x] sequences in group 1, ls[1][x] squences in group 2).
 * \param CL the constraint list
 * \param list_in the diagonals
 * \param n_in number of sequences?
 */
int cl2pair_list_ecl_pc ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in)
{
  int p1, p2, si, s, r, t_s, t_r,t_w, t_s2, t_r2, t_w2, n;
  int a, b, l1, l2;
  int **pos,**list;
  int max_n;

  int nused;
  int *used_list;
  int *sl2, **inv_pos;

  int **nr;


  float nscore, score, tot, filter, avg=0, new=0;
  float **used;

  if ( !A) return 0;
  list=list_in[0];
  n=n_in[0];
  max_n=read_array_size (list, sizeof (int*));

  pos=aln2pos_simple ( A,-1, ns, ls);
  inv_pos=vcalloc ((CL->S)->nseq, sizeof (int*));
  for (a=0; a<ns[1]; a++)inv_pos[ls[1][a]] =seq2inv_pos(A->seq_al[ls[1][a]]);

  l1=strlen (A->seq_al[ls[0][0]]);
  l2=strlen (A->seq_al[ls[1][0]]);
  sl2=vcalloc ((CL->S)->nseq, sizeof (int));
  nr=declare_int (2, MAX(l1,l2)+1);

  for (a=0; a<l1;a++)
    for (b=0; b<ns[0]; b++)
      if (!is_gap(A->seq_al[ls[0][b]][a]))nr[0][a+1]++;
  for (a=0; a<l2;a++)
    for (b=0; b<ns[1]; b++)
      if (!is_gap(A->seq_al[ls[1][b]][a]))nr[1][a+1]++;
  
  for (a=0;a<ns[1]; a++)sl2[ls[1][a]]=1;
  
  CL=index_res_constraint_list (CL, CL->weight_field);
  
  used=declare_float (l2+1,2);
  used_list=vcalloc (l2+1, sizeof (int));
  nused=0;
  
  for (p1=0; p1<=l1; p1++)
    {
      
      for (tot=0,nused=0,si=0;p1>0 && si<ns[0]; si++)
	{
	  s=ls [0][si];r=pos[s][p1-1];
	  for (a=1; r>0 && a<CL->residue_index[s][r][0];a+=3)
	    {
	      t_s=CL->residue_index[s][r][a];t_r=CL->residue_index[s][r][a+1];t_w=CL->residue_index[s][r][a+2];
	      for (b=0; b<CL->residue_index[t_s][t_r][0];)
		{
		  if (b==0){t_s2=t_s;t_r2=t_r;t_w2=t_w;b++;}
		  else 
		    { 
		      t_s2=CL->residue_index[t_s][t_r][b];t_r2=CL->residue_index[t_s][t_r][b+1];t_w2=CL->residue_index[t_s][t_r][b+2];b+=3;
		    }
		  
		  if (sl2[t_s2])
		    {
		      p2=inv_pos[t_s2][t_r2];
		      //score=((float)t_w/(float)NORM_F)*((float)t_w2/(float)NORM_F);
		      score=MIN(((float)t_w/(float)NORM_F),((float)t_w2/(float)NORM_F));
		      
		      if (!used[p2][1] && score>0)
			{
			  used_list[nused++]=p2;
			}
		      
		      tot+=score;
		      used[p2][0]+=score;
		      used[p2][1]++;
		    }
		}
	    }
	}
      //FILTER: Keep in the graph the edges where (p1->p2/(Sum (P1->x))>0.01
      filter=0.01;
      
      for (a=0; a<nused; a++)
	{
	  
	  p2=used_list[a];
	  nscore=used[p2][0]/tot; //Normalized score used for filtering
	  score =used[p2][0];
		  
	  used[p2][0]=used[p2][1]=0;
	  if (n==max_n)
	    {
	      max_n+=10000;list=vrealloc (list, max_n*sizeof (int*));
	    }
	  if (nscore>filter && p1!=0 && p2!=0 && p1!=l1 && p2!=l2)
	    {
	      if (!list[n])
			  list[n]=vcalloc (7, sizeof (int));
	      
	      list[n][0]=p1;
	      list[n][1]=p2;
	      list[n][3]=(l1-(p1))+(p2);
	      score/=(float)((CL->S)->nseq*nr[0][p1]*nr[1][p2]);
	      list[n][2]=(int)((float)score*(float)NORM_F);
	      avg+=(int)((float)score*(float)NORM_F);
	      new++;
	      n++;
	    }
	}
    }
  free_float (used, -1);
  vfree (used_list);
  free_int (inv_pos, -1);
  free_int (pos, -1);
  vfree (sl2);
  free_int (nr, -1);

  n_in[0]=n;
  list_in[0]=list;
  if (new)avg/=new;
  return avg;
}


int cl2pair_list_ecl_rawquad ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in)
{
  int p1, p2, si, s, r, t_s, t_r,t_w, t_s2, t_r2, t_w2,t_s3, t_r3, t_w3, n,n2;
  int a, b, c,l1, l2;
  int score;
  int **pos;
  int **list;
  int max_n;

  int tn;
  
  int nused;
  int *used_list, **used;
  int *sl2, **inv_pos;

  int filter=0;
  int nseq2;
  if ( !A) return 0;
  list=list_in[0];
  n=n_in[0];
  max_n=read_array_size (list, sizeof (int*));
  
 
  n2=0;
  pos=aln2pos_simple ( A,-1, ns, ls);
  inv_pos=vcalloc ((CL->S)->nseq, sizeof (int*));
  for (a=0; a<ns[1]; a++)inv_pos[ls[1][a]] =seq2inv_pos(A->seq_al[ls[1][a]]);
  
  l1=strlen (A->seq_al[ls[0][0]]);
  l2=strlen (A->seq_al[ls[1][0]]);
  sl2=vcalloc ((CL->S)->nseq, sizeof (int));
  
  for (a=0;a<ns[1]; a++)sl2[ls[1][a]]=1;
 
  CL=index_res_constraint_list (CL, CL->weight_field);
  
  used=declare_int (l2+1,2);
  used_list=vcalloc (l2+1, sizeof (int));
  nused=0;
  nseq2=(CL->S)->nseq*(CL->S)->nseq;
  
  for (p1=0; p1<=l1; p1++)
    {
      for (nused=0,si=0;p1>0 && si<ns[0]; si++)
	{
	  s=ls [0][si];r=pos[s][p1-1];
	  for (a=1; r>0 && a<CL->residue_index[s][r][0];a+=3)
	    {
	      t_s=CL->residue_index[s][r][a];t_r=CL->residue_index[s][r][a+1];t_w=CL->residue_index[s][r][a+2];
	      for (b=0; b<CL->residue_index[t_s][t_r][0];)
		{
		  if (b==0){t_s2=t_s;t_r2=t_r;t_w2=t_w;b++;}
		  else 
		    { 
		      t_s2=CL->residue_index[t_s][t_r][b];t_r2=CL->residue_index[t_s][t_r][b+1];t_w2=CL->residue_index[t_s][t_r][b+2];b+=3;
		    }
		  if (sl2[t_s2])
		    {
		      for (c=0; c<CL->residue_index[t_s2][t_r2][0];)
			{
			  if (c==0){t_s3=t_s2;t_r3=t_r2;t_w3=t_w2;c++;}
			  else 
			    { 
			      t_s3=CL->residue_index[t_s2][t_r2][c];t_r3=CL->residue_index[t_s2][t_r2][c+1];t_w3=CL->residue_index[t_s2][t_r2][c+2];c+=3;
			    }
			  
			  if (sl2[t_s3])
			    {
			      p2=inv_pos[t_s3][t_r3];
			      score=MIN(t_w,t_w2);
			      score=MIN(score,t_w3);
			      if (!used[p2][1] && score>0)
				{
				  used_list[nused++]=p2;
				}
			      used[p2][0]+=score;
			      used[p2][1]++;
			    }
			}
		    }
		}
	    }
	}
      for (a=0; a<nused; a++)
	{
	  
	  p2=used_list[a];
	  score=used[p2][0];
	  
	  used[p2][0]=used[p2][1]=0;
	  	  
	  if (n==max_n){max_n+=1000;list=vrealloc (list, max_n*sizeof (int*));}
	  if (score>filter && p1!=0 && p2!=0 && p1!=l1 && p2!=l2)
	    {
	      if (!list[n])list[n]=vcalloc (7, sizeof (int));
	      
	      list[n][0]=p1;
	      list[n][1]=p2;
	      list[n][3]=(l1-(p1))+(p2);
	      list[n][2]=score;
	      n++;
	    }
	}
    }
  

  vfree (used);
  vfree (used_list);
  free_int (inv_pos, -1);
  free_int (pos, -1);
  vfree (sl2);

  n_in[0]=n;
  list_in[0]=list;
  return 1;
}
int cl2pair_list_ecl_norm ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in)
{
  int p1, p2, si, s, r, t_s, t_r, n,n2;
  int a, b, l1, l2;
  int score;
  int **pos;
  int **list;
  int max_n;


  
  int nused;
  int *used_list, *used;
  int *sl2, **inv_pos;

  int filter=0;
  
  if ( !A) return 0;
  list=list_in[0];
  n=n_in[0];
  max_n=read_array_size (list, sizeof (int*));
  
 
  n2=0;
  pos=aln2pos_simple ( A,-1, ns, ls);
  inv_pos=vcalloc ((CL->S)->nseq, sizeof (int*));
  for (a=0; a<ns[1]; a++)inv_pos[ls[1][a]] =seq2inv_pos(A->seq_al[ls[1][a]]);
  
  l1=strlen (A->seq_al[ls[0][0]]);
  l2=strlen (A->seq_al[ls[1][0]]);
  sl2=vcalloc ((CL->S)->nseq, sizeof (int));
  
  for (a=0;a<ns[1]; a++)sl2[ls[1][a]]=1;
 
  CL=index_res_constraint_list (CL, CL->weight_field);
  
  used=vcalloc (l2+1, sizeof (int));
  used_list=vcalloc (l2+1, sizeof (int));
  nused=0;
  
  
  
  for (p1=0; p1<=l1; p1++)
    {
      for (si=0;p1>0 && si<ns[0]; si++)
	{
	  s=ls [0][si];
	  r=pos[s][p1-1];
	  
	  for (a=1; r>0 && a<CL->residue_index[s][r][0];a+=3)
	    {
	      
	      t_s=CL->residue_index[s][r][a];
	      t_r=CL->residue_index[s][r][a+1];
	      
	      for (b=0; b<CL->residue_index[t_s][t_r][0];)
		{
		  int t_s2, t_r2;
		  if (b==0){t_s2=t_s;t_r2=t_r;b++;}
		  else { t_s2=CL->residue_index[t_s][t_r][b];t_r2=CL->residue_index[t_s][t_r][b+1];b+=3;}
		  
		  if (sl2[t_s2])
		    {
		      p2=inv_pos[t_s2][t_r2];
		      if (!used[p2]){used[p2]=1;used_list[nused++]=p2;}
		      else used[p2]++;
		    }
		}
	      
	    }
	}
      if (p1==0 || p1==l1)
	{
	  for (nused=0,p2=0; p2<=l2; p2++)used_list[nused++]=p2;
	}
      else
	{
	  if (!used[0])used_list[nused++]=0;
	  if (!used[l2])used_list[nused++]=l2;
	}
      for (a=0; a<nused; a++)
	{
	  p2=used_list[a];
	  if (p2==0 || p1==0)score=0;
	  
	  else score=(CL->get_dp_cost) (A, pos, ns[0], ls[0], p1-1, pos, ns[1], ls[1],p2-1,CL);	
	  if (score>filter && p1!=0 && p2!=0 && p1!=l1 && p2!=l2)
	    {
	      if (n==max_n){max_n+=1000;list=vrealloc (list, max_n*sizeof (int*));}
	      if (!list[n])list[n]=vcalloc (7, sizeof (int));
	      list[n][0]=p1;
	      list[n][1]=p2;
	      list[n][3]=(l1-(p1))+(p2);
	      list[n][2]=score;
	      n++;
	      if (p1!=0 && p2!=0 && p1!=l1 && p2!=l2)n2++;
	    }
	  used[p2]=0;
	}
      
      nused=0;
    }
  
  vfree (used);
  vfree (used_list);
  free_int (inv_pos, -1);
  free_int (pos, -1);
  vfree (sl2);
  n_in[0]=n;
  list_in[0]=list;
  

   return 1;
}




int cl2pair_list_ref( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list, int *n)
{
  int a, b, l1, l2, n2=0;
  int score;
  int **pos;
  int max_n;
  if (!A) return 0;
  pos=aln2pos_simple ( A,-1, ns, ls);
  l1=strlen (A->seq_al[ls[0][0]]);
  l2=strlen (A->seq_al[ls[1][0]]);
  
  max_n=read_array_size (list[0], sizeof (int));
  
  for (a=0; a<=l1; a++)
    for (b=0; b<=l2; b++)
      {
	score=(a==0 || b==0)?0:slow_get_dp_cost_pc(A, pos, ns[0], ls[0],a-1, pos, ns[1], ls[1], b-1, CL);

	
	if ( score>0 && a!=0 && b!=0 && a!=l1 && b!=l2)
	  {
	    if (n[0]==max_n){max_n+=1000;list[0]=vrealloc (list[0], max_n*sizeof (int*));}
	    if (!list[0][n[0]])list[0][n[0]]=vcalloc (7, sizeof (int));
	    list[0][n[0]][0]=a;
	    list[0][n[0]][1]=b;
	    list[0][n[0]][3]=(l1-a)+b;
	    list[0][n[0]][2]=score;
	    if ( a!=0 && b!=0 && a!=l1 && b!=l2)
	      {
		n2++;
	      }
	    n[0]++;
	  }

      }

  return n[0];
  }

int list2linked_pair_wise ( Alignment *A, int *ns, int **l_s, Constraint_list *CL, int **list, int n);
int two_pass_linked_pair_wise ( Alignment *A, int *ns, int **l_s, Constraint_list *CL)
{
  int n=0, **list=NULL;
  int nscore;
  int mode=2;
  int id;
  
  cl2pair_list (A,ns, l_s, CL, &list, &n, mode, 10);
  nscore=list2linked_pair_wise (A, ns, l_s, CL, list, n);
  id=sub_aln2sim (A, ns, l_s, "idmat_sim");
  
  if (id>50)return nscore;
  ungap_sub_aln ( A, ns[0], l_s[0]);
  ungap_sub_aln ( A, ns[1], l_s[1]);
  cl2pair_list (A,ns, l_s, CL, &list, &n,mode,0);
  nscore=list2linked_pair_wise (A, ns, l_s, CL, list, n);
  cl2pair_list (NULL,ns, l_s, CL, &list, &n, mode, 0);
  return nscore;
}
int clinked_pair_wise ( Alignment *A, int *ns, int **l_s, Constraint_list *CL)
{
  int n=0, **list=NULL;
  int nscore, pscore=0;
  int mode=2;

  int add=0;
  cl2pair_list (A,ns, l_s, CL, &list, &n, mode, 1000);
  nscore=list2linked_pair_wise (A, ns, l_s, CL, list, n);
  HERE ("***********First: %d", nscore);
  if (add)
    {
      while (nscore>pscore)
	{
	  pscore=nscore;
	  ungap_sub_aln ( A, ns[0], l_s[0]);
	  ungap_sub_aln ( A, ns[1], l_s[1]);
	  cl2pair_list (A,ns, l_s, CL, &list, &n, mode, 10);
	  nscore=list2linked_pair_wise (A, ns, l_s, CL, list, n);
	  HERE ("****************New: %d", nscore);
	}
    }
  cl2pair_list (NULL,ns, l_s, CL, &list, &n, mode, 0);
  return nscore;
}
int linked_pair_wise ( Alignment *A, int *nsi, int **lsi, Constraint_list *CL)
{
  int n=0;
  static int **list=NULL;
  int score, a;
  int *ns, **ls;
  int mode=1;//1:ecl, 0:ref

  ns=vcalloc (2, sizeof (int));
  ns[0]=nsi[1]; ns[1]=nsi[0];
  
  ls=declare_int (2, ns[0]+ns[1]);
  for (a=0; a<ns[1]; a++)
    ls[1][a]=lsi[0][a];
  for (a=0; a<ns[0]; a++)
    ls[0][a]=lsi[1][a];
  
  
  cl2pair_list (A,ns, ls, CL, &list, &n, mode, 0);
    
  score=list2linked_pair_wise (A, ns, ls, CL, list, n);
  cl2pair_list (NULL,ns, ls, CL, &list, &n, mode, 0);
  return score;
}
 

int list2linked_pair_wise_nolgp ( Alignment *A, int *ns, int **l_s, Constraint_list *CL, int **list, int n);
int list2linked_pair_wise ( Alignment *A, int *ns, int **l_s, Constraint_list *CL, int **list, int n)
 {
   int mode=0;
   if (mode==0)return list2linked_pair_wise_nolgp    (A, ns, l_s, CL, list, n);
   return 0;
 }


#define LIN(a,b,c) a[b*5+c]
int list2linked_pair_wise_nolgp( Alignment *A, int *ns, int **l_s, Constraint_list *CL, int **list, int n)
{
  int a,b,c, i, j, LEN=0, start_trace;
  int pi, pj,ij, delta_i, delta_j, prev_i, prev_j;
  static int **slist;
  static long *MI, *MJ, *MM,*MT2;
  static int *sortseq;
  static int max_size;
  int gop, gep, igop, igep;
  int l1, l2, l, ls;
  char **al;
  char **aln,*char_buf;
  int ni=0, nj=0;
  long score;
  int nomatch;
 
  l1=strlen (A->seq_al[l_s[0][0]]);
  l2=strlen (A->seq_al[l_s[1][0]]);
  al=declare_char (2,l1+l2+1);
  
 
  //Penalties: max score is NORM_F
  //Penalties must be negative
  igop=CL->gop;
  gep=igep=CL->gep;
  
  if (n>max_size)
    {
      max_size=n;
      
      vfree (MI);vfree (MJ); vfree (MM);
      free_int (slist, -1);
     
      slist=declare_int (n,3);
      
      MI=vcalloc (5*n, sizeof (long));
      MJ=vcalloc (5*n, sizeof (long));
      MM=vcalloc (5*n, sizeof (long));
      
    }
  else
    {
      for (a=0; a<n; a++)
	for (b=0; b<5; b++)LIN(MI,a,b)=LIN(MJ,a,b)=LIN(MJ,a,b)=-1000000;
    }
  
  /*New Bit: Start*/
  if (!sortseq) sortseq=vcalloc( 7, sizeof (int));
  sortseq[0]=0; sortseq[1]=1;sortseq[2]=-1;
  sort_list_int2 (list, sortseq,7, 0, n-1);
  for (a=0; a<n; a++)
  {
      slist[a][0]=a;
      list[a][4]=a;
  }
  
  sortseq[0]=1; sortseq[1]=0;sortseq[2]=-1;
  sort_list_int2 (list, sortseq,7, 0, n-1);
  for (a=0; a<n; a++)
  {
	slist[a][1]=list[a][4];
	list[a][5]=a;
  }
 
  sortseq[0]=3; sortseq[1]=0;sortseq[2]=1;sortseq[3]=-1;
  sort_list_int2 (list, sortseq,7, 0, n-1);
  for (a=0; a<n; a++)
    {
      slist[a][2]=list[a][4];
      list[a][6]=a;
    }

  sortseq[0]=0; sortseq[1]=1;sortseq[2]=-1;
  sort_list_int2 (list, sortseq,7, 0, n-1);
 
  /*New Bit: EnD*/
  
  
 
 
  
  
  for (a=0; a<n; a++)
    {

      
      i=list[a][0];
      j=list[a][1];

      
      if (i==l1 || j==l2)gop=0;
      else gop=igop;

      if (i==l1 && j==l2)start_trace=a;
      else if ( i==0 || j==0)
	{
	  LIN(MM,a,0)=-1000000;
	  if (j==0)
	    {
	      
	      LIN(MJ,a,0)=-10000000;
	      LIN(MI,a,0)=gep*i;
	      
	    }
	  else if (i==0)
	    {
	      
	      LIN(MI,a,0)=-10000000;
	      LIN(MJ,a,0)=gep*j;
	      
	    }

	  LIN(MI,a,1)=LIN(MJ,a,1)=-1;
	  LIN(MI,a,2)=LIN(MJ,a,2)=i;
	  LIN(MI,a,3)=LIN(MJ,a,3)=j;
	  continue;
	}
      
      pi=list[a][5];
      pi=slist[pi-1][1];
      
      pj=list[a][4];
      pj=slist[pj-1][0]; 
      
      ij=list[a][6];
      ij=slist[ij-1][2];
      
      
      ij=list[a][6];
      ij=slist[ij-1][2];
      
      
	
     
      
      prev_i=list[pi][0];
      prev_j=list[pj][1];
      
      delta_i=list[a][0]-list[pi][0];
      delta_j=list[a][1]-list[pj][1];
      
      /*Linear Notation*/
      LIN(MI,a,0)=MAX(LIN(MI,pi,0),(LIN(MM,pi,0)+gop))+delta_i*gep;
      LIN(MI,a,1)=pi;
      LIN(MI,a,2)=delta_i;
      LIN(MI,a,3)=0;
      LIN(MI,a,4)=(LIN(MI,pi,0)>=(LIN(MM,pi,0)+gop))?'i':'m';
    
      
      LIN(MJ,a,0)=MAX(LIN(MJ,pj,0),(LIN(MM,pj,0)+gop))+delta_j*gep;
      LIN(MJ,a,1)=pj;
      LIN(MJ,a,2)=0;
      LIN(MJ,a,3)=delta_j;
      
      LIN(MJ,a,4)=(LIN(MJ,pj,0)>=LIN(MM,pj,0)+gop)?'j':'m';
      
    
      if (a>1 && (ls=list[a][0]-list[ij][0])==(list[a][1]-list[ij][1]))
	{
	  LIN(MM,a,0)=MAX3(LIN(MM,ij,0),LIN(MI,ij,0),LIN(MJ,ij,0))+list[a][2]-(ls*CL->nomatch);

	  LIN(MM,a,1)=ij;
	  LIN(MM,a,2)=ls;
	  LIN(MM,a,3)=ls;
	  if ( LIN(MM,ij,0)>=LIN(MI,ij,0) && LIN(MM,ij,0)>=LIN(MJ,ij,0))LIN(MM,a,4)='m';
	  else if ( LIN(MI,ij,0) >= LIN(MJ,ij,0))LIN(MM,a,4)='i';
	  else LIN(MM,a,4)='j';
	  
	}
      else
	{
	  LIN(MM,a,0)=UNDEFINED;
	  LIN(MM,a,1)=-1;
	}  
    }
  
  a=start_trace;
  if (LIN(MM,a,0)>=LIN(MI,a,0) && LIN(MM,a,0) >=LIN(MJ,a,0))MT2=MM;
  else if ( LIN(MI,a,0)>=LIN(MJ,a,0))MT2=MI;
  else MT2=MJ;

  score=MAX3(LIN(MM,a,0), LIN(MI,a,0), LIN(MJ,a,0));
  
  i=l1;
  j=l2;
  
  
  while (!(i==0 &&j==0))
    {
      int next_a;
      l=MAX(LIN(MT2,a,2),LIN(MT2,a,3));
      // HERE ("%c from %c %d %d SCORE=%d [%d %d] [%2d %2d]", T2[a][5],T2[a][4], T2[a][2], T2[a][3], T2[a][0], gop, gep, i, j);
      if (i==0)
	{
	  while ( j>0)
	    {
	      al[0][LEN]=0;
	      al[1][LEN]=1;
	      j--; LEN++;
	    }
	}
      else if (j==0)
	{
	  while ( i>0)
	    {
	      al[0][LEN]=1;
	      al[1][LEN]=0;
	      i--; LEN++;
	    }
	}
      
      else if (l==0) {HERE ("L=0 i=%d j=%d",l, i, j);exit (0);}
      else 
	{
	  for (b=0; b<l; b++, LEN++)
	    {
	      if (LIN(MT2,a,2)){al[0][LEN]=1;i--;ni++;}
	      else al[0][LEN]=0;
	      
	      if (LIN(MT2,a,3)){al[1][LEN]=1;j--;nj++;}
	      else al[1][LEN]=0;
	    }
	  
	  next_a=LIN(MT2,a,1);
	  if (LIN(MT2,a,4)=='m')MT2=MM;
	  else if (LIN(MT2,a,4)=='i')MT2=MI;
	  else if (LIN(MT2,a,4)=='j')MT2=MJ;
	  a=next_a;
	}
    }
  
 
  
  invert_list_char ( al[0], LEN);
  invert_list_char ( al[1], LEN);
  
  	
  if ( A->declared_len<=LEN)A=realloc_aln  ( A,2*LEN+1);
  aln=A->seq_al;
  char_buf= vcalloc (LEN+1, sizeof (char));
	
  for ( c=0; c< 2; c++)
    {
      for ( a=0; a< ns[c]; a++) 
	{		
	  int ch=0;
	  for ( b=0; b< LEN; b++)
	    {		   
	      if (al[c][b]==1)
		char_buf[b]=aln[l_s[c][a]][ch++];
	      else
		char_buf[b]='-';
	    }
	  char_buf[b]='\0';
	  sprintf (aln[l_s[c][a]],"%s", char_buf);
	}
    }
  
  A->len_aln=LEN;
  A->nseq=ns[0]+ns[1];
  
  vfree (char_buf);
  free_char (al, -1);
  
  return score;
}
int ** aln2local_penalties4link (Alignment *A, int n, int *ls, Constraint_list *CL, int **lg);
int ** aln2local_penalties4link (Alignment *A, int n, int *ls, Constraint_list *CL, int **lg)
{
  //adapted from gap_count in MAFFT V 5.5
  int p,s,l, c1, c2;
  int gep,gop;
  int open=3, close=4, gap=5;
  
  gop=CL->gop;
  gep=CL->gep;
  
  l=strlen (A->seq_al[ls[0]]);
  
  if (!lg)
    {
      lg=declare_int (6, l);
    }
  
  if ( read_array_size_new (lg[0])<l)
    {
      free_int (lg, -1);
      lg=declare_int (6, l);
    }
  
  for( s=0; s<n; s++ ) 
	{
	  c1='x';
	  for (p=0; p<l; p++)
	    {
	      c2=A->seq_al[ls[s]][p];
	      
	      if (c1!='-' && c2=='-')lg[open][p]++;
	      if (c1=='-' && c2!='-')lg[close][p]++;
	      if ( c1=='-')lg[gap][p]++;
	      c1=c2;
	    }
	}
  
  for (p=0; p<l; p++)
    {
      float go, gc, nn;
      nn=n;
      go=lg[open ][p];
      gc=lg[close][p];
     
    
      lg[GOP][p]=0.5*(1-(go/nn))*gop;
      lg[GCP][p]=0.5*(1-(gc/nn))*gop;
      //Checked locacal gep => gives low quality results
      lg[GEP][p]=gep;//(1-((float)lg[gap][p]/(float)n))*gep;
      lg[open][p]=lg[close][p]=lg[gap][p]=0;
      
    }

  return lg;
}

/*********************************COPYRIGHT NOTICE**********************************/
/*© Centro de Regulacio Genomica */
/*and */
/*Cedric Notredame */
/*Tue Oct 27 10:12:26 WEST 2009. */
/*All rights reserved.*/
/*This file is part of T-COFFEE.*/
/**/
/*    T-COFFEE is free software; you can redistribute it and/or modify*/
/*    it under the terms of the GNU General Public License as published by*/
/*    the Free Software Foundation; either version 2 of the License, or*/
/*    (at your option) any later version.*/
/**/
/*    T-COFFEE is distributed in the hope that it will be useful,*/
/*    but WITHOUT ANY WARRANTY; without even the implied warranty of*/
/*    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the*/
/*    GNU General Public License for more details.*/
/**/
/*    You should have received a copy of the GNU General Public License*/
/*    along with Foobar; if not, write to the Free Software*/
/*    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA*/
/*...............................................                                                                                      |*/
/*  If you need some more information*/
/*  cedric.notredame@europe.com*/
/*...............................................                                                                                                                                     |*/
/**/
/**/
/*	*/
/*********************************COPYRIGHT NOTICE**********************************/