#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <stdarg.h>
#include <string.h>
#include <ctype.h> 
#include "io_lib_header.h"
#include "util_lib_header.h"
#include "dp_lib_header.h"
#include "define_header.h"
/************************************************************************************/
/*                NEW      ANALYZE 2    : SAR                                        */
/************************************************************************************/
float display_prediction_old (int **prediction, int n, Alignment *A, Alignment *S, int field);

float display_prediction (int ***count, Alignment *S, int c, int n);
Alignment * filter_aln4sar0 ( Alignment *A, Alignment *S, int c, int leave, char *mode);
Alignment * filter_aln4sar1 ( Alignment *A, Alignment *S, int c, int leave, char *mode);
Alignment * filter_aln4sar2 ( Alignment *A, Alignment *S, int c, int leave, char *mode);
Alignment * filter_aln4sar3 ( Alignment *A, Alignment *S, int c, int leave, char *mode);
Alignment * filter_aln4sar4 ( Alignment *A, Alignment *S, int c, int leave, char *mode);
Alignment * filter_aln4sar5 ( Alignment *A, Alignment *S, int c, int leave, char *mode);

int **sar2profile ( Alignment *A, Alignment *S, int c, int leave);
int **sar2profile_sim ( Alignment *A, Alignment *S, int **sim, int comp, int leave);
int sar_profile2score ( char *seq, int **profile);
double sar_vs_iseq1( char *sar, int  *seq, float gl, int **sim, char *best_aa);
double sar_vs_seq1 ( char *sar, char *seq, float gl, int **sim, char *best_aa);
double sar_vs_seq2 ( char *sar, char *seq, float ng, int **mat, char *a);
double sar_vs_seq3 ( char *sar, char *seq, float ng, int **mat, char *a);

double sar_vs_iseq4 ( char *sar, int *seq, float ng, int **mat, char *a);//supports an extended alphabet
double sar_vs_seq4 ( char *sar, char *seq, float ng, int **mat, char *a);

double sar_vs_seq5 ( char *sar, char *seq, float ng, int **mat, char *a);
int make_sim_pred ( Alignment *A,Alignment *S, int comp, int seq);

int **sar2profile_sim ( Alignment *A, Alignment *S, int **sim, int comp, int leave)
{

  int a, b, r, c, c1, c2, r1, r2, s, p;
  int ***cache, **profile;
  
  
  profile=declare_int (A->len_aln, 26);
  cache=declare_arrayN (3,sizeof (int),2,A->len_aln, 26);
  
  for ( a=0; a< A->len_aln; a++)
    for ( b=0; b< A->nseq; b++)
      {
	r=tolower(A->seq_al[b][a]);
	c=( S->seq_al[comp][b]=='I')?1:0;
	if (b==leave || is_gap(r)) continue;
	cache [c][a][r-'a']++;
      }
  for (a=0; a< A->nseq; a++)
    {
      if ( a==leave) continue;
      for ( b=0; b< A->nseq; b++)
	{
	  c1=(S->seq_al[comp][a]=='I')?1:0;
	  c2=(S->seq_al[comp][b]=='I')?1:0;
	  if ( b==leave || b==a || c1!=1 || c1==c2) continue;
	  s=sim[a][b];
	  
	  for (p=0; p<A->len_aln; p++)
	    {
	      r1=tolower(A->seq_al[a][p]);
	      r2=tolower(A->seq_al[b][p]);
	      if ( is_gap(r1) || is_gap(r2) || r1==r2)continue;
	      r1-='a';r2-='a';
	      if (cache[1][p][r2])continue;
	      if ( s<50)continue;
	      profile[p][r2]-=s;
	    }
	}
    }

  free_arrayN((void***)cache,3);
  return profile;

}
int **sar2profile ( Alignment *A, Alignment *S, int comp, int leave)
{

  int a, b,c,r, n, v, npos=0;
  int ***cache, **profile;
  int ncat;
  float n_gap, max_gap;
  profile=declare_int (A->len_aln, 26);
  cache=declare_arrayN (3,sizeof (int),2,A->len_aln, 26);



  for ( n=0, a=0; a< A->nseq; a++)
    {
      if ( a==leave) continue;
      else n+=(S->seq_al[comp][a]=='I')?1:0;
    }

  for ( a=0; a< A->len_aln; a++)
    for ( b=0; b< A->nseq; b++)
      {
	r=tolower(A->seq_al[b][a]);
	c=( S->seq_al[comp][b]=='I')?1:0;
	if (b==leave) continue;
	else if (is_gap(r))continue;
	r-='a';
	cache [c][a][r]++;
      }

  ncat=15; /*ncat: limit the analysis to columns containing less than ncat categories of aa*/
  max_gap=0.05;
  for (a=0; a< A->len_aln; a++)
    {
      for (n_gap=0,b=0; b< A->nseq; b++)
	n_gap+=(is_gap(A->seq_al[b][a]));
      n_gap/=(float)A->nseq;
      
      if ( n_gap> max_gap)continue;
      
      for (v=0,r=0; r< 26; r++)
	{
	  if (cache [0][a][r] || cache[1][a][r])v++;
	} 
      
      for (n=0,r=0; r< 26 && v<ncat; r++)
	{
	  if (cache [0][a][r] && !cache[1][a][r])
	    {
	      n++;
	      profile[a][r]=-cache[0][a][r];
	    }
	}
      if (n) npos++;
    }
  
  free_arrayN((void***)cache,3);
  return profile;

}
Alignment * filter_aln4sar0 ( Alignment *A, Alignment *S, int comp, int leave, char *mode)
{
  return copy_aln (A,NULL);
}
Alignment * filter_aln4sar1 ( Alignment *inA, Alignment *S, int comp, int leave, char *mode)
{
  Alignment *F, *A;
  int a, b,c, i,r, n0, n1,g, score;
  int ***cache, **list1, **list2;
  int Delta;

  int T1;
  
  /*Keep only the positions where there are residues ONLY associated with 0 sequences*/
  
  list1=declare_int ( inA->nseq, 2);
  list2=declare_int ( inA->len_aln, 2);

  cache=declare_arrayN (3,sizeof (int),inA->len_aln,2, 26);
  F=copy_aln (inA, NULL);
  
  A=copy_aln (inA, NULL);
  A->nseq=strlen (S->seq_al[comp]);

  strget_param (mode, "_T1_", "5", "%d", &T1);
  for ( a=0; a< A->len_aln; a++)
    {
      n1=n0=g=0;
      for (b=0; b< A->nseq; b++)
	{
	  if ( b==leave) continue;
	  i=(S->seq_al[comp][b]=='I')?1:0;
	  r=tolower(A->seq_al[b][a]);
	  if ( r=='-')continue;
	  cache[a][i][r-'a']++;
	}
    }
  
  for (a=0; a< A->nseq; a++)
    for ( score=0,b=0; b<A->len_aln; b++)
      {
	r=tolower (A->seq_al[a][b]);
	if ( is_gap(r))continue;
	else if ( cache[b][0][r-'a'] && !cache[b][1][r-'a'])list1[a][0]++;
      }
 
  for (a=0; a< A->len_aln; a++)
    {
      for ( score=0,b=0; b< A->nseq; b++)
	{
	  r=tolower (A->seq_al[b][a]);
	  if ( r=='-')continue;
	  else r-='a';
	  if ( cache[a][0][r] && !cache[a][1][r])score ++;
	}
      list2[a][0]=a;
      list2[a][1]=score;
    }
  sort_int (list2, 2, 1, 0, F->len_aln-1);
  
  Delta=A->len_aln/(100/T1);
  for ( a=0; a< F->len_aln-Delta; a++)
    {
      b=list2[a][0];
      for ( c=0; c<F->nseq; c++)
	{
	  F->seq_al[c][b]='-';
	}
    }

  ungap_aln (F);
  free_aln (A);
  free_arrayN ( (void ***)cache, 3);
  free_arrayN ((void**)list1, 2);
  free_arrayN ((void**)list2, 2);

  return F;
}
Alignment * filter_aln4sar2 ( Alignment *inA, Alignment *S, int comp, int leave, char *mode)
{
  Alignment *F, *A;
  int a,b,r,ncat;
  int *cache;
  int max_ncat=10;

  /*Keep Low entropy columns that contain less than ncat categories of different amino acids*/
  /*REmove columns containing 10% or more gaps*/
  
  cache=vcalloc ( 500, sizeof (char));
  F=copy_aln (inA, NULL);
  A=copy_aln (inA, NULL);
  A->nseq=strlen (S->seq_al[comp]);
  for ( a=0; a< A->len_aln; a++)
    {
      for (ncat=0,b=0; b< A->nseq; b++)
	{
	  if ( b==leave) continue;

	  r=tolower(A->seq_al[b][a]);
	  if ( !cache[r])ncat++;
	  cache[r]++;
	}
      
      if ( ncat <max_ncat && ((cache['-']*100)/A->nseq)<10)
	{
	  ;
	}
      else
	{
	  for (b=0; b<F->nseq; b++)
	    {
	      r=tolower(F->seq_al[b][a]);
	      F->seq_al[b][a]='-';
	      cache[r]=0;
	    }
	}
      for (b=0; b<A->nseq; b++)
	  {
	    r=tolower(A->seq_al[b][a]);
	    cache[r]=0;
	  }
    }

  free_aln (A);
  ungap_aln (F);
  vfree (cache);
  return F;
}

Alignment * filter_aln4sar3 ( Alignment *inA, Alignment *S, int comp, int leave, char *mode)
{
  Alignment *F, *rA, *A;
  int a, b,c;
  int **list1;
  char *bufS, *bufA;
  int Delta;
  int T3;
  
  /*Keep the 10% positions most correlated with the 0/1 pattern*/
  
  A=copy_aln (inA, NULL);
  A->nseq=strlen (S->seq_al[comp]);
  F=copy_aln (inA, NULL);
  rA=rotate_aln (A, NULL);
  
  strget_param (mode, "_T3_", "10", "%d", &T3);
  
  
  list1=declare_int ( inA->len_aln, 2);
  bufA=vcalloc ( A->nseq+1, sizeof (char));
  bufS=vcalloc ( A->nseq+1, sizeof (char));
  
  sprintf ( bufS, "%s", S->seq_al[comp]);
  splice_out_seg(bufS,leave, 1);
  
  
  for (a=0; a< A->len_aln; a++)
    {
      char aa;
      list1[a][0]=a;
      sprintf (bufA, "%s", rA->seq_al[a]);
      splice_out_seg (bufA,leave,1);
      list1[a][1]=(int)sar_vs_seq3 ( bufS, bufA,0,NULL, &aa);
    }

  sort_int (list1, 2, 1, 0, F->len_aln-1);
  Delta=F->len_aln/(100/T3);
  for ( a=0; a< F->len_aln-Delta; a++)
    {
	  b=list1[a][0];
	  
	  for ( c=0; c<F->nseq; c++)
	    {
	      F->seq_al[c][b]='-';
	    }

    }
  F->score_aln=list1[F->len_aln-1][1];
  ungap_aln (F);

  free_aln (rA);    
  free_aln(A);
  free_arrayN ((void**)list1, 2);
  vfree (bufS);vfree (bufA);
  return F;
}
Alignment * filter_aln4sar4 ( Alignment *inA, Alignment *S, int comp, int leave, char *mode)
{
  Alignment *F, *A;
  int a, b,c, i,r, n0, n1,g,score;
  int ***cache, **list1, **list2;
 
  /*Keep only the positions where there are residues ONLY associated with 0 sequences*/
  
  list1=declare_int ( inA->nseq, 2);
  list2=declare_int ( inA->len_aln, 2);

  cache=declare_arrayN (3,sizeof (int),inA->len_aln,2, 26);
  F=copy_aln (inA, NULL);
  A=copy_aln (inA, NULL);
  A->nseq=strlen (S->seq_al[comp]);
  
  for ( a=0; a< A->len_aln; a++)
    {
      n1=n0=g=0;
      for (b=0; b< A->nseq; b++)
	{
	  if ( b==leave) continue;
	  i=(S->seq_al[comp][b]=='I')?1:0;
	  r=tolower(A->seq_al[b][a]);
	  if ( r=='-')continue;
	  cache[a][i][r-'a']++;
	  n1+=i;
	}
    }
  
  
  for (a=0; a< A->len_aln; a++)
    {
      for ( score=0,b=0; b< A->nseq; b++)
	{
	  r=tolower (F->seq_al[b][a]);
	  if ( r=='-')continue;
	  else r-='a';
	  if (cache[a][1][r]>=n1/2)score=1;
	}
      list2[a][0]=a;
      list2[a][1]=score;
    }

  
  for ( a=0; a< F->len_aln; a++)
    {
      if ( list2[a][1]==1);
      else
	{
	  b=list2[a][0];
	  for ( c=0; c<F->nseq; c++)
	    {
	      F->seq_al[c][b]='-';
	    }
	}
    }
  ungap_aln (F);
  free_aln (A);
  free_arrayN ( (void ***)cache, 3);
  free_arrayN ((void**)list1, 2);
  free_arrayN ((void**)list2, 2);

  return F;
}

Alignment * filter_aln4sar5 ( Alignment *inA, Alignment *S, int comp, int leave, char *mode)
{
  Alignment *F, *rA, *A;
  int a, b,c;
  int **list1;
  char *bufS, *bufA;
  int max;
  /*Look for the positions that show the best correlation between the sequence variation and the SAR*/

  A=copy_aln (inA, NULL);
  A->nseq=strlen (S->seq_al[comp]);
  
  rA=rotate_aln (inA, NULL);
  F=copy_aln (inA, NULL);
  
  list1=declare_int ( A->len_aln, 2);
  bufA=vcalloc ( A->nseq+1, sizeof (char));
  bufS=vcalloc ( A->nseq+1, sizeof (char));
  


  sprintf ( bufS, "%s", S->seq_al[comp]);
  splice_out_seg(bufS,leave, 1);
  
  
  for (a=0; a< A->len_aln; a++)
    {
      char aa;
      list1[a][0]=a;
      sprintf (bufA, "%s", rA->seq_al[a]);
      splice_out_seg (bufA,leave,1);
      list1[a][1]=(int)sar_vs_seq4 ( bufS, bufA,0,NULL, &aa);
    }

  sort_int (list1, 2, 1, 0, F->len_aln-1);
  max=F->score=list1[F->len_aln-1][1];
  max-=(max/10);
  
  
  for ( a=0; a< F->len_aln-10; a++)
    {
      
	  b=list1[a][0];
	  
	  for ( c=0; c<F->nseq; c++)
	    {
	      F->seq_al[c][b]='-';
	    }

    }
  F->score_aln=10;
  ungap_aln (F);
  free_aln (inA);
  free_aln (rA);    
  free_arrayN ((void**)list1, 2);
  vfree (bufS);vfree (bufA);
  return F;
}

int sar_profile2score ( char *seq, int **P)
{
  int a,r, l, score;
  
  l=strlen (seq);
  for ( score=0,a=0; a< l; a++)
    {
      r=seq[a];
      if ( is_gap(r))continue;
      score+=P[a][tolower(r)-'a'];
    }
  return score;
}
int make_sim_pred ( Alignment *A,Alignment *S, int comp, int seq)
{
  int a, b, i, r1, r2;
  static float **cscore;
  static float **tscore;

  if ( !cscore)
    {
      cscore=declare_float (2, 2);
      tscore=declare_float (2, 2);
    }
  
  for (a=0; a< 2; a++)for (b=0; b<2; b++)cscore[a][b]=tscore[a][b]=0;
  
  for ( a=0; a<A->len_aln; a++)
    {
      r1=A->seq_al[seq][a];
      if ( r1=='-') continue;
      else
	{
	  for ( b=0; b< A->nseq; b++)
	    {
	      if (b==seq) continue;
	      else
		{
		  r2=A->seq_al[b][a];
		  if (r2=='-')continue;
		  else
		    {
		      
		      i=(S->seq_al[comp][b]=='I')?1:0;
		      cscore[i][0]+=(r1==r2)?1:0;
		      cscore[i][1]++;
		    }
		}
	    }
	 
	  for (i=0; i<2; i++)
	    {
	      cscore[i][0]/=(cscore[i][1]==0)?1:cscore[i][1];
	      tscore[i][0]+=cscore[i][0];tscore[i][1]++;
	      cscore[i][0]=cscore[i][1]=0;
	    }
	}
    }

  fprintf ( stdout, "\nn\t 1: %.2f 0: %.2f", tscore[1][0],tscore[0][0]);
  return ( tscore[1][0]>=tscore[0][0])?1:0;
}
      

Alignment * sar_analyze (Alignment *inA, Alignment *inS, char *mode)
{
  int ***sim,***glob_results, ***comp_results;
  int *count;
  int a,b,c,m;
  float *tot2;
  Alignment *A=NULL,*S=NULL,*F, *SUBSET;
  char *subset, *target;
  int jack, T, filter;
  filter_func *ff;
  int n_methods=0;
  char *prediction, *reliability;
  int pred_start=0, pred_end, ref_start=0, ref_end;
  int display, CSV=1, NONCSV=0;
  char method[5];
  
  strget_param (mode, "_METHOD_", "1111", "%s_", method);
  ff=vcalloc (6,sizeof (filter_func));
  if (method[0]=='1')ff[n_methods++]=filter_aln4sar0;
  if (method[1]=='1')ff[n_methods++]=filter_aln4sar1;
  if (method[2]=='1')ff[n_methods++]=filter_aln4sar2;
  if (method[3]=='1')ff[n_methods++]=filter_aln4sar3;
  /*
    ff[n_methods++]=filter_aln4sar4;
    ff[n_methods++]=filter_aln4sar5;
  */
  sim=vcalloc (n_methods, sizeof (int**));
  

  tot2=vcalloc ( 10, sizeof (float));
  subset=vcalloc ( 100, sizeof (char));
  target=vcalloc ( 100, sizeof (char));
  
  strget_param (mode, "_TARGET_", "no", "%s_", target);
  strget_param (mode, "_SUBSET_", "no", "%s_", subset);
  strget_param (mode, "_JACK_", "0", "%d", &jack);
  strget_param (mode, "_T_", "0", "%d", &T);
  strget_param (mode, "_FILTER_", "11", "%d", &filter);
  strget_param (mode, "_DISPLAY_", "0", "%d", &display);
  
  

  if ( !strm (target, "no"))
    {
      Alignment *T;
      T=main_read_aln(target, NULL);
      if ( T->len_aln !=inA->len_aln )
	{
	  printf_exit ( EXIT_FAILURE,stderr, "Error: %s is incompatible with the reference alignment [FATAL:%s]",target,PROGRAM);
	}
      
      inA=stack_aln (inA, T);
      
    }

  if ( !strm(subset, "no")) 
    {
      SUBSET=main_read_aln (subset, NULL);
      sarset2subsarset ( inA, inS, &A, &S, SUBSET);
    }
  else
    {
      A=inA;
      S=inS;
    }
  

  prediction=vcalloc ( n_methods+1, sizeof (char));
  reliability=vcalloc ( n_methods+1, sizeof (char));
  
  glob_results=declare_arrayN(3, sizeof (int), n_methods*2, 2, 2);

  count=vcalloc (S->nseq, sizeof (int));
  for (a=0; a<S->nseq; a++)
    {
      int l;
      l=strlen (S->seq_al[a]);
      for ( b=0; b<l; b++)
	count[a]+=(S->seq_al[a][b]=='I')?1:0;
    }
  if ( display==CSV)
    {fprintf ( stdout, "\nCompound %s ; Ntargets %d", S->name[a],count[a]);
      pred_start=(strlen (S->seq_al[0])==A->nseq)?0:strlen (S->seq_al[0]);
      pred_end=A->nseq;
      for (a=pred_start; a< pred_end; a++)
	fprintf ( stdout, ";%s", A->name[a]);
      fprintf ( stdout, ";npred;");
    }
  
  
  for (a=0; a<S->nseq; a++)
    {
      int n_pred;
      comp_results=declare_arrayN(3, sizeof (int), n_methods*2, 2, 2);
      
      pred_start=(strlen (S->seq_al[a])==A->nseq)?0:strlen (S->seq_al[a]);
      pred_end=A->nseq;
      if ( display==CSV)fprintf ( stdout, "\n%s;%d", S->name[a],count[a]);
      
      for (n_pred=0,b=pred_start; b<pred_end;b++)
	{
	  int t, score=0,pred, real;
	  
	  if ( display==NONCSV)fprintf ( stdout, "\n>%-15s %10s %c ", S->name[a], A->name[b], (pred_start==0)?S->seq_al[a][b]:'?');
	  if (jack || b==pred_start)
	    {
	      for (m=0; m<n_methods; m++)
		{
		  free_int (sim[m], -1);
		  F=(ff[m]) (A,S,a,(jack==0)?-1:b, mode);
		  sim[m]=aln2sim_mat(F, "idmat");
		  free_aln (F);
		}
	    }
	  
	  for (m=0; m<n_methods; m++)
	    {
	      int Nbsim=0,Ybsim=0,bsim=0;
	      ref_start=0;
	      ref_end=strlen (S->seq_al[m]);
	      
	      for (c=ref_start;c<ref_end; c++)
		{
		  if ( b==c) continue;
		  else if ( S->seq_al[a][c]=='O')
		    {
		      Nbsim=MAX(Nbsim,sim[m][b][c]);
		    }
		  else 
		    {
		      Ybsim=MAX(Ybsim,sim[m][b][c]);
		    }
		}
	      
	      bsim=(Ybsim>Nbsim)?Ybsim:-Nbsim;
	      pred=(bsim>0)?1:0;
	      real=(S->seq_al[a][b]=='O')?0:1;
	      comp_results[m][pred][real]++;
	      glob_results[m][pred][real]++;
	      score+=pred;
	      prediction[m]=pred+'0';
	      reliability[m]=(FABS((Ybsim-Nbsim))-1)/10+'0';
	    }
	  
	  if ( score>0)n_pred++;
	  prediction[m]=reliability[m]='\0';
	  if (display==NONCSV)fprintf ( stdout, "Compound_Count:%d primary_predictions: %s Total: %d", count[a],prediction, score);
	  else if ( display==CSV)fprintf ( stdout, ";%d", score);
	  for (t=0; t<n_methods; t++)
	    {
	      if (score>t)
		{
		  comp_results[t+n_methods][1][real]++;
		  glob_results[t+n_methods][1][real]++;
		}
	      else 
		{
		  comp_results[t+n_methods][0][real]++;
		  glob_results[t+n_methods][0][real]++;
		}
	    }
	}
      if ( display==NONCSV)
	{if ( pred_start==0)display_prediction (comp_results, S,a, n_methods*2);}
      else fprintf (stdout, ";%d;",n_pred);
    }
  if ( display==NONCSV)if (pred_start==0)display_prediction (glob_results, S,-1, n_methods*2);
      
  
  exit (EXIT_SUCCESS);
}
float display_prediction (int ***count, Alignment *S, int c, int n)
{
  float tp,tn,fn,fp,sp,sn,sn2;
  int a, nm;

  nm=n/2;
  
  for (a=0; a<n; a++)
    {
      tp=count[a][1][1];
      tn=count[a][0][0];
      fp=count[a][1][0];
      fn=count[a][0][1];

      sn2=tp/(tp+fp);
      sn=tp/(tp+fn);
      sp=tn/(tn+fp);
      if ( a<nm)fprintf ( stdout, "\n>#Method %d Compound %15s sp=%.2f sn=%.2f sn2=%.2f",a, (c==-1)?"TOTAL":S->name[c],sp, sn, sn2 );
      else fprintf ( stdout, "\n>#Combined: T=%d Compound %15s sp=%.2f sn=%.2f sn2=%.2f",a-nm, (c==-1)?"TOTAL":S->name[c],sp, sn, sn2 );
    }
  fprintf ( stdout, "\n");
  return 0; 
}
 
float display_prediction_2 (int **prediction, int n,Alignment *A, Alignment *S, int field)
{
  int a, t, T;
  float max_sn, max_sp;
  
  if ( field==17 || field ==18) 
    {
      printf_exit ( EXIT_FAILURE, stderr, "\nERROR: Do not use filed %d in display_prediction", field);
    }
  
  sort_int_inv ( prediction, 10,field, 0, n-1);
  for (t=0,a=0; a<n; a++)
    {
      t+=prediction[a][3];
      prediction[a][17]=t;
    }

  for (t=0,a=n-1; a>=0; a--)
    {
      prediction[a][18]=t;
      t+=prediction[a][3];
    }

  max_sn=max_sp=T=0;
  for (a=0; a<n; a++)
    {
      float tp, fn, fp, sp, sn;
      
      tp=prediction[a][17];
      fn=prediction[a][18];
      fp=(a+1)-tp;
      
      sp=((tp+fp)==0)?0:tp/(tp+fp);
      sn=((tp+fn)==0)?0:tp/(tp+fn);

      if (sp>0.8)
	{
	  if (sn>max_sn)
	    {
	      max_sn=sn;
	      max_sp=sp;
	      
	      T=prediction[a][field];
	    }
	}
    }
  if (max_sn==0)
      fprintf (stdout, "\n T =%d SN=%.2f SP= %.2f",T,max_sn,max_sp);
  else
      fprintf (stdout, "\n T =%d SN=%.2f SP= %.2f",T,max_sn,max_sp);
  
  return max_sn;
}


/************************************************************************************/
/*                NEW      ANALYZE     : SAR                                        */
/************************************************************************************/
float** cache2pred1 (Alignment *A,int**cache, int *ns, int **ls, Alignment *S, char *compound, char *mode);
float** cache2pred2 (Alignment *A,int**cache, int *ns, int **ls, Alignment *S, char *compound, char *mode);
float** cache2pred3 (Alignment *A,int**cache, int *ns, int **ls, Alignment *S, char *compound, char *mode);
float** cache2pred4 (Alignment *A,int**cache, int *ns, int **ls, Alignment *S, char *compound, char *mode);
float** cache2pred5 (Alignment *A,int**cache, int *ns, int **ls, Alignment *S, char *compound, char *mode);
float** cache2pred_new (Alignment *A,int**cache, int *ns, int **ls, Alignment *S, char *compound, char *mode);

int **sar2cache_adriana ( Alignment *A, int *ns,int **ls, Alignment *S, char *compound, char *mode);
int **sar2cache_proba_old ( Alignment *A, int *ns,int **ls, Alignment *S, char *compound, char *mode);
int **sar2cache_count1 ( Alignment *A, int *ns,int **ls, Alignment *S, char *compound, char *mode);
int **sar2cache_count2 ( Alignment *A, int *ns,int **ls, Alignment *S, char *compound, char *mode);
int **sar2cache_count3 ( Alignment *A, int *ns,int **ls, Alignment *S, char *compound, char *mode);

int **sar2cache_proba_new ( Alignment *A, int *ns,int **ls, Alignment *S, char *compound, char *mode);
int **sar2cache_proba2 ( Alignment *A, int *ns,int **ls, Alignment *S, char *compound, char *mode);
int **analyze_sar_compound1 ( char *name, char *seq, Alignment *A, char *mode);
int **analyze_sar_compound2 ( char *name, char *seq, Alignment *A, char *mode);

int aln2n_comp_col ( Alignment *A, Alignment *S, int ci);




int ***simple_sar_analyze_vot ( Alignment *inA, Alignment *SAR, char *mode);
int ***simple_sar_analyze_col ( Alignment *inA, Alignment *SAR, char *mode);



int sarset2subsarset ( Alignment *A, Alignment *S, Alignment **subA, Alignment **subS, Alignment *SUB);
int benchmark_sar (int v);
int aln2jack_group1 (Alignment *A, int seq, int **l1, int *nl1, int **l2, int *nl2);
int aln2jack_group2 (Alignment *A, int seq, int **l1, int *nl1, int **l2, int *nl2);
int aln2jack_group3 (Alignment *A, char *sar_seq, int **l1, int *nl1, int **l2, int *nl2);
float** jacknife5 (Alignment*A,int **cache, int *ns, int **ls, Alignment *S, char *compound, char *mode);
float** jacknife6 (Alignment*A,int **cache, int *ns, int **ls, Alignment *S, char *compound, char *mode);

int process_cache ( Alignment *A,Alignment *S, int ***Cache, char *mode);
Alignment *analyze_compounds (Alignment *A, Alignment *S, char *mode);

Alignment *analyze_compounds (Alignment *A, Alignment *S, char *mode)
{
  int a, b, c, tot, n;
  int **sim;
  int sar1, sar2;
  
  sim=aln2sim_mat (A, "idmat");
  for (a=0; a< S->nseq; a++)
    {
      for (n=0, tot=0, b=0; b< A->nseq-1; b++)
	{
	  sar1=(S->seq_al[a][b]=='I')?1:0;
	  for ( c=b+1; c<A->nseq; c++)
	    {
	      sar2=(S->seq_al[a][c]=='I')?1:0;
	      
	      if (sar1 && sar2)
		{
		  tot+=sim[b][c];
		  n++;
		}
	    }
	}
      fprintf ( stdout, ">%-10s   CMPSIM: %.2f\n", S->name[a],(float)tot/(float)n); 
    }
  free_int (sim, -1);
  return A;
}

int print_seq_pos ( int pos, Alignment *A, char *seq);
int abl1_evaluation (int p);
int print_seq_pos ( int pos, Alignment *A, char *seq)
{
  int a, b, s;
  
  s=name_is_in_list (seq, A->name, A->nseq, MAXNAMES);
  fprintf ( stdout, "S=%d", s);
  
  for (b=0,a=0; a<pos; a++)
    {
      if (!is_gap (A->seq_al[s][a]))b++;
    }
  fprintf ( stdout, "Pos %d SEQ %s: %d ", pos+1, seq, b+246);
  if ( strm ( seq, "ABL1")) fprintf ( stdout , "PT: %d", abl1_evaluation (b+246));
  return 0;
}

int process_cache ( Alignment *A,Alignment *S, int  ***Cache, char *mode)
{
  int a, b;
  int **pos, **pos2;
  int **C;
  int ab1, *ab1_pos;
  int weight_mode;
  
  strget_param ( mode, "_WEIGHT_", "1", "%d", &weight_mode);
  pos=declare_int(A->len_aln+1,2);
  pos2=declare_int (A->len_aln+1,S->nseq);
  for (a=0; a<S->nseq; a++)
    {
      C=Cache[a];
      for (b=0; b< A->len_aln; b++)
	{
	    pos[b][0]+=C[26][b];
	    if ( C[26][b]>0)
	      {
		pos[b][1]++;
		pos2[b][a]=1;
	      }
	}
    }
  
  C=Cache[0];
  ab1=name_is_in_list ("ABL1", A->name, A->nseq,100);
  ab1_pos=vcalloc (A->len_aln+1, sizeof (int));
  
  for ( b=0,a=0; a< A->len_aln; a++)
    {
      if ( A->seq_al[ab1][a]=='-')ab1_pos[a]=-1;
      else ab1_pos[a]=++b;
    }
    
  for ( a=0; a< A->len_aln; a++)
    {
      fprintf ( stdout, "\n%4d %5d %5d %5d [%c] [%2d] ALN", a+1, pos[a][0], pos[a][1], ab1_pos[a]+246,A->seq_al[ab1][a],abl1_evaluation (ab1_pos[a]+246));
      for ( b=0; b< S->nseq; b++)fprintf ( stdout, "%d", pos2[a][b]);
    }
  return 1;
}
int abl1_evaluation (int p)
{
  if ( p==248) return 10;
  if ( p==250) return 10;
  if ( p==253) return 10;
  if ( p==254) return 10;
  if ( p==255) return 9;
  if ( p==256) return 10;
  if ( p==257) return 5;
  if ( p==258) return 8;
  if ( p==269) return 8;
  if ( p==291) return 4;
  if ( p==294) return 8;
  if ( p==299) return 10;
  if ( p==306) return 0;
  if ( p==314) return 9;
  if ( p==315) return 10;
  if ( p==318) return 10;

  if ( p==319) return 10;
  if ( p==321) return 10;
  if ( p==323) return 0;
  if ( p==324) return 0;
  if ( p==339) return 0;
  if ( p==340) return 0;
  if ( p==355) return 5;
  if ( p==364) return 10;
  
  if ( p==366) return 0;
  if ( p==368) return 10;
  if ( p==370) return 10;
  if ( p==372) return 0;
  if ( p==378) return 8;
  if ( p==382) return 10;

  if ( p==384) return 10;
  if ( p==387) return 10;
  if ( p==395) return 8;

  if ( p==398) return 8;
  if ( p==399) return 8;
  if ( p==400) return 8;
  if ( p==403) return 0;
  if ( p==416) return 8;
  if ( p==419) return 5;
  if ( p>400) return 0;
  return -1;
}
float** cache2pred1 (Alignment*A,int **cache, int *ns, int **ls, Alignment *S, char *compound, char *mode)
{
  int s1, s2, seq1, seq2, r1, r2,col, pred, real, ci;
  double score, max, id, m;
  float **R, T;
  
  
  int used_col, used_res,is_used_col, n_res=0;
  int weight_mode;
  /*Predict on ns[1] what was trained on ns[0]*/

  strget_param ( mode, "_THR_", "0.09", "%f", &T);
  strget_param ( mode, "_WEIGHT_", "0", "%d", &weight_mode);
  
  R=declare_float (2, 2);
  ci=name_is_in_list ( compound, S->name, S->nseq, -1);

  

  for (s1=0; s1<ns[1]; s1++)
    {
      int v;
      seq1=ls[1][s1];
      
      for (max=0,score=0, col=0; col<A->len_aln; col++)
	{
	  int max1;
	  r1=tolower (A->seq_al[seq1][col]);
	  for (max1=0,id=0, m=0,s2=0; s2<ns[0]; s2++)
	    {
	      seq2=ls[0][s2];
	      if ( S->seq_al[ci][seq2]=='O')continue;
	      if ( cache[seq2][col]==0 && !is_gap( A->seq_al[seq2][col]))continue;	      
	      
	      r2=tolower ( A->seq_al[seq2][col]);
	      if ( is_gap(r2))continue;
	      
	      v=(cache[seq2][col]>0 && weight_mode==1)?cache[seq2][col]:1;

	      max+=v;
	      if ( r2==r1)
		{
		  score+=v;
		}
	      
	    }
	  
	}
      pred=(( score/max) >T)?1:0;
      real=(S->seq_al[ci][seq1]=='I')?1:0;
      R[pred][real]++;
      
      fprintf ( stdout, "\n>%s %d%d SCORE %.2f C %s [SEQ]\n", A->name[seq1],real, pred, (float)score/(float)max, compound);
    }

  for (used_col=0,used_res=0,col=0; col<A->len_aln; col++)
    {
      for (is_used_col=0,s2=0; s2<ns[0]; s2++)
	{
	  seq2=ls[0][s2];
	  if ( cache[seq2][col]==0 && !is_gap(A->seq_al[seq2][col]))n_res++;
	  else if (is_gap(A->seq_al[seq2][col]));
	  else 
	    {
	    is_used_col=1;
	    used_res++;
	    }
	}
      used_col+=is_used_col;
    }
  fprintf ( stdout, "\n>%s USED_POSITIONS: COL: %.2f RES: %.2f COMP\n", S->name[ci],  (float)used_col/(float)A->len_aln, (float)used_res/(float) n_res);
  
  return R;
}

float** cache2pred2 (Alignment*A,int **cache, int *ns, int **ls, Alignment *S, char *compound, char *mode)
{
  int s1, s2, seq1, seq2, r1, r2,col, pred, real, ci;
  double score, max;
  float **R, T;
  
  
  int used_col, used_res,is_used_col, n_res=0;
  /*Predict on ns[1] what was trained on ns[0]*/

  strget_param ( mode, "_THR_", "0.5", "%f", &T);
  
  
  R=declare_float (2, 2);
  ci=name_is_in_list ( compound, S->name, S->nseq, -1);

  for (s1=0; s1<ns[1]; s1++)
    {
      int v;
      seq1=ls[1][s1];
      fprintf ( stdout, "\n");
      for (max=0,score=0, col=0; col<A->len_aln; col++)
	{
	  int used;
	  
	  r1=tolower (A->seq_al[seq1][col]);
	  for (used=0,s2=0; s2<ns[0]; s2++)
	    {
	      seq2=ls[0][s2];
	      
	      if ( S->seq_al[ci][seq2]=='O')continue;
	      if ( cache[seq2][col]==0 && !is_gap( A->seq_al[seq2][col]))continue;	      
	      

	      r2=tolower ( A->seq_al[seq2][col]);
	      if ( is_gap(r2))continue;
	      
	      v=cache[seq2][col];
	      if ( r2==r1){score+=v;}
	      used=1;
	      max+=v;
	    }
	  if (used) fprintf ( stdout, "%c", r1);
	}

      pred=(( score/max) >T)?1:0;
      real=(S->seq_al[ci][seq1]=='I')?1:0;
      R[pred][real]++;

      fprintf ( stdout, "PSEQ: %-10s SC: %4d MAX: %4d S: %.2f R: %4d", A->name[seq1],(int)score, (int)max, (float)score/max,real);

    }

  for (used_col=0,used_res=0,col=0; col<A->len_aln; col++)
    {
      for (is_used_col=0,s2=0; s2<ns[0]; s2++)
	{
	  seq2=ls[0][s2];
	  if ( cache[seq2][col]==0 && !is_gap(A->seq_al[seq2][col]))n_res++;
	  else if (is_gap(A->seq_al[seq2][col]));
	  else 
	    {
	    is_used_col=1;
	    used_res++;
	    }
	}
      used_col+=is_used_col;
    }
  fprintf ( stdout, "\n>%s USED_POSITIONS: COL: %.2f RES: %.2f COMP\n", S->name[ci],  (float)used_col/(float)A->len_aln, (float)used_res/(float) n_res);
  
  return R;
}

float** cache2pred3 (Alignment*A,int **cache, int *ns, int **ls, Alignment *S, char *compound, char *mode)
{
  int s1, s2, seq1, seq2, r1, r2,col, pred, real, ci, a, n;
  double score, max;
  float **R, T;
  
  
  
  int tp, tn, fn, fp;
  int best_tp, best_fp;
  int delta, best_delta;
  int **list;
  
  /*Predict on ns[1] what was trained on ns[0]*/

  strget_param ( mode, "_THR_", "0.5", "%f", &T);
  
  
  R=declare_float (2, 2);
  ci=name_is_in_list ( compound, S->name, S->nseq, -1);
  list=declare_int ( ns[1],3);
  
  for (s1=0; s1<ns[1]; s1++)
    {
      int v;
      seq1=ls[1][s1];
      
      for (max=0,score=0, col=0; col<A->len_aln; col++)
	{
	  int used;
	  
	  r1=tolower (A->seq_al[seq1][col]);
	  for (used=0,s2=0; s2<ns[0]; s2++)
	    {
	      seq2=ls[0][s2];
	      
	      if ( S->seq_al[ci][seq2]=='O')continue;
	      if ( cache[seq2][col]==0 && !is_gap( A->seq_al[seq2][col]))continue;	      
	      

	      r2=tolower ( A->seq_al[seq2][col]);
	      if ( is_gap(r2))continue;
	      
	      v=cache[seq2][col];
	      if ( r2==r1){score+=v;}
	      used=1;
	      max+=v;
	    }
	}

     

      pred=(( score/max) >T)?1:0;
      real=(S->seq_al[ci][seq1]=='I')?1:0;
      
      list[s1][0]=real;
      list[s1][1]=(int)((score/max)*(float)1000);
      list[s1][2]=seq1;
      
      

    }
  sort_int_inv (list, 3, 1, 0, ns[1]-1);
    
  for ( a=0; a<ns[1]; a++)
    {
      seq1=list[a][2];
      fprintf ( stdout, "PSEQ: %-10s SC: %5d R: %4d\n", A->name[seq1],list[a][0], list[a][1]);
    }

  for (n=0, a=0; a<ns[1]; a++)n+=list[a][0];
  for (best_delta=100000,best_tp=0,tp=0,fp=0,best_fp=0,a=0; a<ns[1]; a++)
    {
      tp+=list[a][0];
      fp+=1-list[a][0];
      delta=(n-(tp+fp));
      if (FABS(delta)<best_delta)
	{
	  best_delta=delta;
	  best_tp=tp;
	  best_fp=fp;
	}
    }
  /*R[pred][real]*/
  tp=best_tp;
  fp=best_fp;
  fn=n-tp;
  tn=ns[1]-(tp+fp+fn);
  R[1][1]=tp;
  R[1][0]=fp;
  R[0][1]=fn;
  R[0][0]=tn;
  free_int (list, -1);
  return R;
}
float** cache2pred4 (Alignment*A,int **cache, int *ns, int **ls, Alignment *S, char *compound, char *mode)
{
  int s1, s2, seq1, seq2, ci, a,b, c, n;
  double score;
  float **R;


  int tp, tn, fn, fp;
  int best_tp, best_fp;
  int delta, best_delta;
  int **list;
  int **sim;
  int *ul;
  int nused=0;
  
  /*Predict on ns[1] what was trained on ns[0]*/
  /*Identify interesting coloumns*/
  ul=vcalloc ( A->len_aln, sizeof (int));
  for (a=0; a< A->len_aln; a++)
    for ( b=0; b< A->nseq; b++)
      if ( cache[b][a])ul[nused++]=a;
  
  /*compute the similarity on the used columns*/
  
  R=declare_float (2, 2);
  sim=declare_int (A->nseq, A->nseq);
  for (a=0; a< A->nseq; a++)
    for ( b=0; b< A->nseq; b++)
      {
	for (c=0; c< nused; c++)
	  {
	    if ( A->seq_al[a][ul[c]]==A->seq_al[b][ul[c]])sim[a][b]++;
	  }
	sim[a][b]=(sim[a][b]*100)/nused;
      }
  vfree (ul);
  
    
  
  
  ci=name_is_in_list ( compound, S->name, S->nseq, -1);
  list=declare_int ( ns[1],2);
  
  for (s1=0; s1<ns[1]; s1++)
    {

      seq1=ls[1][s1];

      for (score=0,s2=0; s2<ns[0]; s2++)
	{
	  seq2=ls[0][s2];

	  if ( seq1==seq2)continue;
	  if (S->seq_al[ci][seq2]=='I')score=MAX(score, sim[seq1][seq2]);
	}
      list[s1][0]=(S->seq_al[ci][seq1]=='I')?1:0;
      list[s1][1]=(int)score;

    }
  sort_int_inv (list, 2, 1, 0, ns[1]-1);
  
  for (n=0, a=0; a<ns[1]; a++)n+=list[a][0];
  for (best_delta=100000,best_tp=0,tp=0,fp=0,best_fp=0,a=0; a<ns[1]; a++)
    {
      tp+=list[a][0];
      fp+=1-list[a][0];
      delta=(n-(tp+fp));
      if (FABS(delta)<best_delta)
	{
	  best_delta=delta;
	  best_tp=tp;
	  best_fp=fp;
	}
    }
  /*R[pred][real]*/
  tp=best_tp;
  fp=best_fp;
  fn=n-tp;
  tn=ns[1]-(tp+fp+fn);
  R[1][1]=tp;
  R[1][0]=fp;
  R[0][1]=fn;
  R[0][0]=tn;
  free_int (list, -1);
  free_int (sim, -1);
  return R;
}

float** cache2pred5 (Alignment*A,int **cache, int *ns, int **ls, Alignment *S, char *compound, char *mode)
{
  int s1, s2, seq1, seq2, ci, a, n;
  double score;
  float **R;
  
  
  
  int tp, tn, fn, fp;
  int best_tp, best_fp;
  int delta, best_delta;
  int **list;
  static int **sim;
  
  /*Predict on ns[1] what was trained on ns[0]*/
  
  R=declare_float (2, 2);

  if ( sim==NULL)
    sim=aln2sim_mat (A, "idmat");
  
  
  
  ci=name_is_in_list ( compound, S->name, S->nseq, -1);
  list=declare_int ( ns[1],2);
  
  for (s1=0; s1<ns[1]; s1++)
    {

      seq1=ls[1][s1];

      for (score=0,s2=0; s2<ns[0]; s2++)
	{
	  seq2=ls[0][s2];

	  if ( seq1==seq2)continue;
	  if (S->seq_al[ci][seq2]=='I')score=MAX(score, sim[seq1][seq2]);
	}
      list[s1][0]=(S->seq_al[ci][seq1]=='I')?1:0;
      list[s1][1]=(int)score;

    }
  sort_int_inv (list, 2, 1, 0, ns[1]-1);
  
  for (n=0, a=0; a<ns[1]; a++)n+=list[a][0];
  for (best_delta=100000,best_tp=0,tp=0,fp=0,best_fp=0,a=0; a<ns[1]; a++)
    {
      tp+=list[a][0];
      fp+=1-list[a][0];
      delta=(n-(tp+fp));
      if (FABS(delta)<best_delta)
	{
	  best_delta=delta;
	  best_tp=tp;
	  best_fp=fp;
	}
    }
  /*R[pred][real]*/
  tp=best_tp;
  fp=best_fp;
  fn=n-tp;
  tn=ns[1]-(tp+fp+fn);
  R[1][1]=tp;
  R[1][0]=fp;
  R[0][1]=fn;
  R[0][0]=tn;
  free_int (list, -1);
  return R;
}

float** jacknife5 (Alignment*A,int **cacheIN, int *ns, int **ls, Alignment *S, char *compound, char *mode)
{
  int seq1, ci, a,b, c, n;
  double score, max_score;
  float **R;


  int tp, tn, fn, fp;
  int best_tp, best_fp;
  int delta, best_delta;
  int **list;
  int **cache;
 
  ci=name_is_in_list ( compound, S->name, S->nseq, -1);
  list=declare_int (A->nseq,2);
  R=declare_float (2, 2);
  

  for ( a=0; a<A->nseq; a++)
    {
      int real, res;

      ns[0]=A->nseq-1;
      ns[1]=1;
      for (c=0,b=0; b<A->nseq; b++)
	if (a!=b)ls[0][c++]=b;
      ls[1][0]=a;
      
      
      cache=sar2cache_count1 (A, ns, ls,S, compound, mode);
      for (b=0; b<=26; b++)
	for ( c=0; c< A->len_aln; c++)
	  cacheIN[b][c]+=cache[b][c];
      
      seq1=a;
      real=(S->seq_al[ci][seq1]=='I')?1:0;
      fprintf ( stdout, ">%-10s %d ", A->name[seq1], real);
      

      
      for (max_score=0,b=0; b<A->len_aln; b++)
	max_score+=cache[26][b];
      
      for (score=0,b=0; b<A->len_aln; b++)
	{
	  res=tolower (A->seq_al[seq1][b]);
	  if ( cache[26][b]==0) continue;
	  if ( !is_gap(res))
	    {
	      score+=cache[res-'a'][b];
	    }
	  /*fprintf ( stdout, "%c[%3d]", res,b);*/
	}
      fprintf ( stdout, " SCORE: %5d SPRED %d RATIO: %.2f \n", (int)score, a, (score*100)/max_score);
      list[a][0]=real;
      
      if ( strstr (mode, "SIMTEST"))list[a][1]=(score*100)/max_score;
      else list[a][1]=(score*100)/max_score;
      free_int (cache, -1);
    }
  

  sort_int_inv (list, 2, 1, 0, A->nseq-1);
  for (n=0, a=0; a<A->nseq; a++)
    {
      n+=list[a][0];
    }

  for (best_delta=100000,best_tp=0,tp=0,fp=0,best_fp=0,a=0; a<A->nseq; a++)
    {
      
      tp+=list[a][0];
      fp+=1-list[a][0];
      delta=(n-(tp+fp));
      if (FABS(delta)<best_delta)
	{
	  best_delta=delta;
	  best_tp=tp;
	  best_fp=fp;
	}
    }
  
  /*R[pred][real]*/
  tp=best_tp;
  fp=best_fp;
  fn=n-tp;
  tn=A->nseq-(tp+fp+fn);
  R[1][1]=tp;
  R[1][0]=fp;
  R[0][1]=fn;
  R[0][0]=tn;
  free_int (list, -1);
  
  return R;
}
float** jacknife6 (Alignment*A,int **cache, int *ns, int **ls, Alignment *S, char *compound, char *mode)
{
  int seq1, ci, a,b, c,d,e,f, n;
  double score;
  float **R;
  
  
  int tp, tn, fn, fp;
  int best_tp, best_fp;
  int delta, best_delta;
  int **list;
  
  ci=name_is_in_list ( compound, S->name, S->nseq, -1);
  list=declare_int (A->len_aln,2);
  R=declare_float (2, 2);
  

  for ( a=0; a<A->nseq; a++)
    {
      int sar, res;
      int **new_cache;
      
      ns[0]=A->nseq-1;
      ns[1]=1;
      for (c=0,b=0; b<A->nseq; b++)
	if (a!=b)ls[0][c++]=b;
      ls[1][0]=a;
      
      cache=sar2cache_proba_new (A, ns, ls,S, compound, mode);
      

      new_cache=declare_int (27,A->len_aln);
      
      for (d=0; d< A->len_aln; d++)
	{
	  int **analyze;
	  if ( cache[26][d]==0)continue;
	  analyze=declare_int (26, 2);
	  
	  for ( e=0; e< ns[0]; e++)
	    {
	      f=ls[0][e];
	      sar=(S->seq_al[ci][f]=='I')?1:0;
	      res=tolower (A->seq_al[f][d]);
	      
	      if ( res=='-') continue;
	      analyze[res-'a'][sar]++;
	    }
	  for (e=0;e<26; e++)
	    {
	      if ( analyze[e][1]){new_cache[26][d]=1;new_cache[e][d]+=cache[e][d];}
	      /*
	      if ( analyze[e][0] && analyze[e][1]){new_cache[26][d]=1;new_cache[e][d]+=analyze[e][1];}
	      else if ( analyze[e][0]){new_cache[26][d]=1;new_cache[e][d]-=analyze[e][0]*10;}
	      else if ( analyze[e][1]){new_cache[26][d]=1;new_cache[e][d]+=analyze[e][1];}
	      else if ( !analyze[e][0] &&!analyze[e][1]);
	      */
	    }
	  free_int (analyze, -1);
	}

      seq1=a;
      sar=(S->seq_al[ci][seq1]=='I')?1:0;
      fprintf ( stdout, ">%-10s %d ", A->name[seq1], sar);
      
      for (score=0,b=0; b<A->len_aln; b++)
	{
	  res=tolower (A->seq_al[seq1][b]);
	  if ( cache[26][b]==0) continue;
	  if ( !is_gap(res))
	    {
	      score+=new_cache[res-'a'][b];
	    }
	}
      fprintf ( stdout, " SCORE: %5d SPRED\n", (int)score);
      list[seq1][0]=sar;
      list[seq1][1]=(int)score;
    
      free_int (new_cache, -1);
      free_int (cache, -1);
    }
  sort_int_inv (list, 2, 1, 0, A->nseq-1);
  for (n=0, a=0; a<A->nseq; a++)n+=list[a][0];
  for (best_delta=100000,best_tp=0,tp=0,fp=0,best_fp=0,a=0; a<A->nseq; a++)
    {
      tp+=list[a][0];
      fp+=1-list[a][0];
      delta=(n-(tp+fp));
      if (FABS(delta)<best_delta)
	{
	  best_delta=delta;
	  best_tp=tp;
	  best_fp=fp;
	}
    }
  
  fprintf ( stderr, "\n%d %d", best_tp, best_fp);
  /*R[pred][real]*/
  tp=best_tp;
  fp=best_fp;
  fn=n-tp;
  tn=A->nseq-(tp+fp+fn);
  R[1][1]=tp;
  R[1][0]=fp;
  R[0][1]=fn;
  R[0][0]=tn;
  free_int (list, -1);
  

  return R;
}
float** cache2pred_new (Alignment*A,int **cache, int *ns, int **ls, Alignment *S, char *compound, char *mode)
{
  int s1, seq1, ci, a,b, n;
  double score;
  float **R;


  int tp, tn, fn, fp;
  int best_tp, best_fp;
  int delta, best_delta;
  int **list;
  
  ci=name_is_in_list ( compound, S->name, S->nseq, -1);
  list=declare_int ( ns[1],2);
  R=declare_float (2, 2);
  
  for (s1=0; s1<ns[1]; s1++)
    {
      int res, real;
     
      seq1=ls[1][s1];
      real=(S->seq_al[ci][seq1]=='I')?1:0;
      fprintf ( stdout, ">%-10s %d ", A->name[seq1], real);
      for (score=0,b=0; b<A->len_aln; b++)
	{
	  res=tolower (A->seq_al[seq1][b]);
	  if ( cache[26][b]==0) continue;
	  if ( !is_gap(res))
	    {
	      score+=cache[res-'a'][b];
	    }
	  fprintf ( stdout, "%c", res);
	}
      fprintf ( stdout, " SCORE: %5d SPRED\n", (int)score);
      list[s1][0]=real;
      list[s1][1]=(int)score;
    }
    
  sort_int_inv (list, 2, 1, 0, ns[1]-1);
  
  for (n=0, a=0; a<ns[1]; a++)n+=list[a][0];
  for (best_delta=100000,best_tp=0,tp=0,fp=0,best_fp=0,a=0; a<ns[1]; a++)
    {
      tp+=list[a][0];
      fp+=1-list[a][0];
      delta=(n-(tp+fp));
      if (FABS(delta)<best_delta)
	{
	  best_delta=delta;
	  best_tp=tp;
	  best_fp=fp;
	}
    }
  
  

  /*R[pred][real]*/
  tp=best_tp;
  fp=best_fp;
  fn=n-tp;
  tn=ns[1]-(tp+fp+fn);
  R[1][1]=tp;
  R[1][0]=fp;
  R[0][1]=fn;
  R[0][0]=tn;
  free_int (list, -1);
  

  return R;
}
float** cache2pred_forbiden_res (Alignment*A,int **cache, int *ns, int **ls, Alignment *S, char *compound, char *mode)
{
  int s1,seq1, ci, a,b, c, n;
  double score;
  float **R;
  
  
  int tp, tn, fn, fp;
  int best_tp, best_fp;
  int delta, best_delta;
  int **list;
  int **new_cache;
  int **mat;

  mat=read_matrice ( "blosum62mt");
  
  ci=name_is_in_list ( compound, S->name, S->nseq, -1);
  list=declare_int ( ns[1],2);
  R=declare_float (2, 2);
  
  for (s1=0; s1<ns[1]; s1++)
    {
      int res, real;
     
      seq1=ls[1][s1];
      real=(S->seq_al[ci][seq1]=='I')?1:0;
      fprintf ( stdout, ">%-10s %d ", A->name[seq1], real);
      for (score=0,b=0; b<A->len_aln; b++)
	{
	  res=tolower (A->seq_al[seq1][b]);
	  if ( cache[26][b]==0) continue;
	  if ( !is_gap(res))
	    {
	      score+=cache[res-'a'][b];
	    }
	  fprintf ( stdout, "%c", res);
	}
      fprintf ( stdout, " SCORE: %5d SPRED\n", (int)score);
      list[s1][0]=real;
      list[s1][1]=(int)score;
    }
  new_cache=declare_int (27,A->len_aln);
  for (a=0; a< A->len_aln; a++)
    {
      int **analyze, real, res, d;
      int *res_type;
      int **sub;
      int *keep;
      keep=vcalloc ( 26, sizeof (int));
      res_type=vcalloc ( 26, sizeof (int));
      sub=declare_int (256, 2);
      
      if ( cache[26][a]==0)continue;
      analyze=declare_int (26, 2);
      for ( b=0; b< ns[0]; b++)
	{
	  seq1=ls[0][b];
	  real=(S->seq_al[ci][seq1]=='I')?1:0;
	  res=tolower (A->seq_al[seq1][a]);
	  
	  if ( res=='-') continue;
	  analyze[res-'a'][real]++;
	}
      fprintf ( stdout, "RSPRED: ");
      for (c=0;c<26; c++)fprintf ( stdout, "%c", c+'a');
      fprintf ( stdout, "\nRSPRED: ");
      for (c=0;c<26; c++)
	{
	  if ( analyze[c][0] && analyze[c][1]){fprintf ( stdout, "1");res_type[c]='1';}
	  else if ( analyze[c][0]){new_cache[26][a]=1;new_cache[c][a]-=analyze[c][0];fprintf ( stdout, "0");res_type[c]='0';}
	  else if ( analyze[c][1]){new_cache[26][a]=1;new_cache[c][a]+=analyze[c][1];fprintf ( stdout, "1");res_type[c]='1';}
	  else if ( !analyze[c][0] &&!analyze[c][1]){fprintf ( stdout, "-");res_type[c]='-';}
	}

     
      for ( c=0; c<26; c++)
	{
	  for ( d=0; d<26; d++)
	    {
	      
	      if ( res_type[c]==res_type[d])
		{
		  sub[res_type[c]][0]+=mat[c][d];
		  sub[res_type[c]][1]++;
		}
	      if ( res_type[c]!='-' && res_type[d]!='-')
		{
		  sub['m'][0]+=mat[c][d];
		  sub['m'][1]++;
		}
	    }
	}
      for ( c=0; c< 256; c++)
	{
	  if ( sub[c][1])fprintf ( stdout, " %c: %5.2f ", c, (float)sub[c][0]/(float)sub[c][1]);
	}
      fprintf ( stdout, " SC: %d\nRSPRED  ", cache[26][a]);
      
      for ( c=0; c<26; c++)
	if ( res_type[c]=='1')
	  {
	    for (d=0; d<26; d++)
	      if (mat[c][d]>0)keep[d]++;
	    keep[c]=9;
	  }

      for (c=0; c<26; c++)
	{
	  if ( keep[c]>10)fprintf ( stdout, "9");
	  else fprintf ( stdout, "%d", keep[c]);
	}
      for ( c=0; c<26; c++)
	{
	  if ( keep[c]>8)new_cache[c][a]=10;
	  else new_cache[c][a]=-10;
	}
      fprintf ( stdout, "\n");
      free_int (analyze, -1);
      free_int (sub, -1);
      vfree (res_type);
      vfree (keep);
      
    }
  for ( a=0; a<25; a++)
    for (b=a+1; b<26; b++)
      {
	int r1, r2;
	r1=a+'a';r2=b+'a';
	if ( strchr("bjoxz", r1))continue;
	if ( strchr("bjoxz",r2))continue;
	
	if ( mat[a][b]>0 && a!=b)fprintf ( stdout, "\nMATANALYZE %c %c %d", a+'a', b+'a', mat[a][b]);
      }
  
  for (s1=0; s1<ns[1]; s1++)
    {
      int res, real;
     
      seq1=ls[1][s1];
      real=(S->seq_al[ci][seq1]=='I')?1:0;
      fprintf ( stdout, ">%-10s %d ", A->name[seq1], real);
      for (score=0,b=0; b<A->len_aln; b++)
	{
	  res=tolower (A->seq_al[seq1][b]);
	  if ( cache[26][b]==0) continue;
	  if ( !is_gap(res))
	    {
	      score+=new_cache[res-'a'][b];
	    }
	  fprintf ( stdout, "%c", res);
	}
      fprintf ( stdout, " SCORE: %5d SPRED\n", (int)score);
      list[s1][0]=real;
      list[s1][1]=(int)score;
    }
  free_int (new_cache, -1);
  sort_int_inv (list, 2, 1, 0, ns[1]-1);
  
  
  for (n=0, a=0; a<ns[1]; a++)n+=list[a][0];
  for (best_delta=100000,best_tp=0,tp=0,fp=0,best_fp=0,a=0; a<ns[1]; a++)
    {
      tp+=list[a][0];
      fp+=1-list[a][0];
      delta=(n-(tp+fp));
      if (FABS(delta)<best_delta)
	{
	  best_delta=delta;
	  best_tp=tp;
	  best_fp=fp;
	}
    }



  /*R[pred][real]*/
  tp=best_tp;
  fp=best_fp;
  fn=n-tp;
  tn=ns[1]-(tp+fp+fn);
  R[1][1]=tp;
  R[1][0]=fp;
  R[0][1]=fn;
  R[0][0]=tn;
  free_int (list, -1);
  

  return R;
}

int **sar2cache_proba_old ( Alignment *A, int *ns,int **ls, Alignment *S, char *compound, char *mode)
{
  int col, s, seq,ms,mseq, res, mres, res1, n,maxn1, maxn2,maxn3, t, ci, a;
  float quant=0;
  int **list;

  int N1msa,N1sar, N, N11, N10, N01,N00, SCORE, COL_INDEX, RES;
  int nfield=0;
  int value;
  float T1, T2, T3, T4;
  int weight_mode;
  int **cache;
  static int **sim;
  int sim_weight, w, sw_thr;
  int train_mode;
  
  float zscore;
  
  RES=nfield++;COL_INDEX=nfield++;N1msa=nfield++;N1sar=nfield++;N=nfield++;N11=nfield++;N10=nfield++;N01=nfield++;N00=nfield++;SCORE=nfield++;
  ci=name_is_in_list ( compound, S->name, S->nseq, -1);
  cache=declare_int (A->nseq, A->len_aln);
  
  strget_param ( mode, "_FILTER1_", "0"   , "%f", &T1);
  strget_param ( mode, "_FILTER2_", "1000000", "%f", &T2);
  strget_param ( mode, "_FILTER3_", "0"   , "%f", &T3);
  strget_param ( mode, "_FILTER4_", "1000000", "%f", &T4);
  strget_param ( mode, "_SIMWEIGHT_", "1", "%d", &sim_weight);
  strget_param ( mode, "_SWTHR_", "30", "%d", &sw_thr);
  strget_param (mode, "_TRAIN_","1", "%d", &train_mode);
  strget_param (mode, "_ZSCORE_","0", "%f", &zscore);
  
  



  if (sim_weight==1 && !sim) sim=aln2sim_mat(A, "idmat");
  for ( ms=0; ms<ns[0]; ms++)
    {
      mseq=ls[0][ms];
      if ( S->seq_al[ci][mseq]!='I')continue;

      list=declare_int (A->len_aln+1, nfield);
      for (t=0,n=0, col=0; col< A->len_aln; col++)
	{
	  int same_res;
	  
	  mres=tolower(A->seq_al[mseq][col]);
	  list[col][RES]=mres;
	  list[col][COL_INDEX]=col;

	  if ( is_gap(mres))continue;
	  for ( s=0; s<ns[0]; s++)
	    {
	      seq=ls[0][s];
	      res=tolower(A->seq_al[seq][col]);
	      if (is_gap(res))continue;
	      

	      if (sim_weight==1)
		{
		  w=sim[seq][mseq];w=(mres==res)?100-w:w;
		  if (w<sw_thr)w=0;
		}
	      else
		w=1;
	      
	      if ( train_mode==4)
		{
		  if ( S->seq_al[ci][seq]=='I')same_res=1;
		  else same_res=(res==mres)?1:0;
		}
	      else
		same_res=(res==mres)?1:0;
	      
	      list[col][N]+=w;
	      
	      if (S->seq_al[ci][seq]=='I' && same_res)list[col][N11]+=w;
	      else if (S->seq_al[ci][seq]=='I' && same_res)list[col][N10]+=w;
	      else if (S->seq_al[ci][seq]=='O' && same_res)list[col][N01]+=w;
	      else if (S->seq_al[ci][seq]=='O' && same_res)list[col][N00]+=w;
	      
	      if ( S->seq_al[ci][seq]=='I')list[col][N1sar]+=w;
	      if ( same_res)list[col][N1msa]+=w;
	      
	    }
	  
	  list[col][SCORE]=(int)evaluate_sar_score1 (list[col][N], list[col][N11], list[col][N1msa], list[col][N1sar]);
	  
	}

      strget_param ( mode, "_MAXN1_", "5", "%d", &maxn1);
      strget_param ( mode, "_WEIGHT_", "1", "%d", &weight_mode);
      strget_param ( mode, "_QUANT_", "0.0", "%f", &quant);
      
      sort_int_inv (list,nfield,SCORE,0,A->len_aln-1);
      if ( quant !=0)
	{
	
	  n=quantile_rank ( list,SCORE, A->len_aln,quant);
	  sort_int (list,nfield,N1msa, 0, n-1);
	  maxn1=MIN(n,maxn1);
	}
      
      for (a=0; a<maxn1; a++)
	{
	  col=list[a][COL_INDEX];
	  res1=list[a][RES];
	  value=list[a][SCORE];
	  if ( value>T1 && value<T2){cache[mseq][col]= value;}
	}
      free_int (list, -1);
    }
  
  /*Filter Columns*/
  list=declare_int (A->len_aln+1, nfield);
  for ( col=0; col< A->len_aln; col++)
    {
      list[col][COL_INDEX]=col;
      for ( s=0; s<ns[0]; s++)
	{
	  seq=ls[0][s];
	  list[col][SCORE]+=cache[seq][col];
	}
    }
 
  /*Filter Columns with a score not between T2 and T3*/
  
  for (col=0; col< A->len_aln; col++)
    if (list[col][SCORE]<T3 || list[col][SCORE]>T4)
      {
	list[col][SCORE]=0;
	for (s=0; s< A->nseq; s++)
	  if (!is_gap(A->seq_al[s][col]))cache[s][col]=0;
      }
  
  /*Keep The N Best Columns*/
  if ( zscore!=0)
    {
      double sum=0, sum2=0, z;
      int n=0;
      for (a=0; a< A->len_aln; a++)
	{
	  if ( list[a][SCORE]>0)
	    {
	      sum+=list[a][SCORE];
	      sum2+=list[a][SCORE]*list[a][SCORE];
	      n++;
	    }
	}
      for (a=0; a<A->len_aln; a++)
	{
	  if ( list[a][SCORE]>0)
	    {
	      z=return_z_score (list[a][SCORE], sum, sum2,n);
	      if ((float)z<zscore)
		{  
		  col=list[a][COL_INDEX];
		  for (s=0; s<A->nseq; s++)
		    cache [s][col]=0;
		}
	      else
		{
		  fprintf ( stdout, "\nZSCORE: KEEP COL %d SCORE: %f SCORE: %d\n", list[a][COL_INDEX], (float)z, list[a][SCORE]);
		}
	    }
	}
    }
  else
    {
      sort_int_inv (list,nfield,SCORE,0,A->len_aln-1);
      strget_param ( mode, "_MAXN2_", "100000", "%d", &maxn2);
      
      for (a=maxn2;a<A->len_aln; a++)
	{
	  col=list[a][COL_INDEX];
	  for (s=0; s<A->nseq; s++)
	    cache [s][col]=0;
	}
    }

  /*Get Rid of the N best Columns*/;
  strget_param ( mode, "_MAXN3_", "0", "%d", &maxn3);
  
  for (a=0; a<maxn3;a++)
    {
      col=list[a][COL_INDEX];
      for (s=0; s<A->nseq; s++)
	cache [s][col]=0;
    }
  
  return cache;
}
int aln2n_comp_col ( Alignment *A, Alignment *S, int ci)
{
  int  res, seq,sar, col, r;
  int **analyze;

  int tot=0;
  
  analyze=declare_int (27, 2);  
  for ( col=0; col< A->len_aln; col++)
    {
      int n1, n0;

      
      for ( n1=0, n0=0,seq=0; seq<A->nseq; seq++)
	{
	  res=tolower(A->seq_al[seq][col]);
	  sar=(S->seq_al[ci][seq]=='I')?1:0;
	  n1+=(sar==1)?1:0;
	  n0+=(sar==0)?1:0;
	  if ( res=='-')continue;
	  res-='a';
	  analyze[res][sar]++;
	}
      
      for (r=0; r<26; r++)
	{
	  int a0,a1;
	  a0=analyze[r][0];
	  a1=analyze[r][1];
	  
	  
	  if ( a1==n1 && a0<n0)
	    {
	      tot++;
	    }
	}
      for ( r=0; r<26; r++)analyze[r][0]=analyze[r][1]=0;
    }
  
  free_int (analyze, -1);
  return tot;
}
int **sar2cache_count1 ( Alignment *A, int *ns,int **ls, Alignment *S, char *compound, char *mode)
{
  int  maxn2, res, seq,sar, ci, col,s, r;
  int **analyze, **list, **cache;
  static int **mat;

  int a0,a1, w;
  if (!mat) mat=read_matrice ("blosum62mt");
  
  
  list=declare_int ( A->len_aln, 2);
  cache=declare_int ( 27, A->len_aln);
  analyze=declare_int (27, 2);  
  
  ci=name_is_in_list ( compound, S->name, S->nseq, -1);
  
  for ( col=0; col< A->len_aln; col++)
    {
      int n1, n0;

      
      for ( n1=0, n0=0,s=0; s<ns[0]; s++)
	{
	  seq=ls[0][s];
	  res=tolower(A->seq_al[seq][col]);
	  sar=(S->seq_al[ci][seq]=='I')?1:0;
	  n1+=(sar==1)?1:0;
	  n0+=(sar==0)?1:0;
	  if ( res=='-')continue;
	  res-='a';
	  	  
	  analyze[res][sar]++;
	}
      
      for (r=0; r<26; r++)
	{
	  
	  a0=analyze[r][0];
	  a1=analyze[r][1];
	  
	  if ( strstr (mode, "SIMTEST"))
	    {
	      w=a1;
	    }
	  else if (a1 )
	    {
	      w=n0-a0;
	    }
	  else w=0;
	  
	  cache[r][col]+=w;
	  cache[26][col]=MAX(w, cache[26][col]);
	}
      
      for ( r=0; r<26; r++)analyze[r][0]=analyze[r][1]=0;
      list[col][0]=col;
      list[col][1]=cache[26][col];
    }
  
  free_int (analyze, -1);

  sort_int_inv (list, 2, 1, 0, A->len_aln-1);
  
  strget_param ( mode, "_MAXN2_", "100000", "%d", &maxn2);

  for ( col=maxn2; col<A->len_aln; col++)
    for ( r=0; r<=26; r++)cache[r][list[col][0]]=0;
 
  free_int (list, -1);
  return cache;
}


int **sar2cache_count2 ( Alignment *A, int *ns,int **ls, Alignment *S, char *compound, char *mode)
{
  int maxn2, res, seq,sar, ci, col,s, r;
  int **analyze, **list, **cache, **conseq;
  static int **mat;
  int w=0;
  if (!mat) mat=read_matrice ("blosum62mt");
  
  
  list=declare_int ( A->len_aln, 2);
  cache=declare_int ( 27, A->len_aln);
  conseq=declare_int ( A->len_aln,3);

  analyze=declare_int (27, 2);  
  
  ci=name_is_in_list ( compound, S->name, S->nseq, -1);  
  for ( col=0; col< A->len_aln; col++)
    {
      int n1, n0;

      for ( n1=0, n0=0,s=0; s<ns[0]; s++)
	{
	  seq=ls[0][s];
	  res=tolower(A->seq_al[seq][col]);
	  sar=(S->seq_al[ci][seq]=='I')?1:0;
	  n1+=(sar==1)?1:0;
	  n0+=(sar==0)?1:0;
	  if ( res=='-')continue;
	  res-='a';
	  analyze[res][sar]++;
	}
      for (r=0; r<26; r++)
	{
	  int a0,a1;
	  a0=analyze[r][0];
	  a1=analyze[r][1];
	  if ( a1==n1 && a0<n0)
	    {
	      	      	      
	      w=n0-a0;
	      conseq[col][0]=r;
	      conseq[col][1]=w;
	    }
	}
      for ( r=0; r<26; r++)analyze[r][0]=analyze[r][1]=0;
    }
  free_int (analyze, -1);
  
  for (s=0; s<ns[0]; s++)
    {
      int w1, w2;
      seq=ls[0][s];
      for (w1=0,w2=0,col=0; col<A->len_aln; col++)
	{

	  res=tolower(A->seq_al[seq][col]);
	  if ( is_gap(res))continue;
	  else res-='a';
	  
	  if ( conseq[col][1] && res!=conseq[col][0])w1++;
	  if ( conseq[col][1])w2++;
	}
      for (col=0; col<A->len_aln; col++)
	{
	  res=tolower(A->seq_al[seq][col]);
	  if ( is_gap(res))continue;
	  else res-='a';
	  
	  if ( conseq[col][1] && res!=conseq[col][0])conseq[col][2]+=(w2-w1);
	}
    }
  
  for (col=0; col<A->len_aln; col++)
    {
      r=conseq[col][0];
      w=conseq[col][2];

      
      cache[r][col]=cache[26][col]=list[col][1]=w;
      list[col][0]=col;
    }
  sort_int_inv (list, 2, 1, 0, A->len_aln-1);
  strget_param ( mode, "_MAXN2_", "100000", "%d", &maxn2);

  for ( col=maxn2; col<A->len_aln; col++)
    for ( r=0; r<=26; r++)cache[r][list[col][0]]=0;
 

  free_int (list, -1);
  return cache;
}  

int **sar2cache_count3 ( Alignment *A, int *ns,int **ls, Alignment *S, char *compound, char *mode)
{
  int  maxn2, res, seq,sar, ci, col,s, r, a1, a0, n1, n0;
  int **analyze, **list, **cache;
  static int **mat;
  
  if (!mat) mat=read_matrice ("blosum62mt");
  
  
  list=declare_int ( A->len_aln, 2);
  cache=declare_int ( 27, A->len_aln);
  analyze=declare_int (27, 2);  
  
  ci=name_is_in_list ( compound, S->name, S->nseq, -1);
  
  for ( col=0; col< A->len_aln; col++)
    {
      double e, g;
      for ( n1=0, n0=0,s=0; s<ns[0]; s++)
	{
	  seq=ls[0][s];
	  res=tolower(A->seq_al[seq][col]);
	  sar=(S->seq_al[ci][seq]=='I')?1:0;
	  n1+=(sar==1)?1:0;
	  n0+=(sar==0)?1:0;
	  if ( res=='-')continue;
	  res-='a';
	  	  
	  analyze[res][sar]++;
	}
      
      /*Gap*/
      for (g=0,r=0; r<A->nseq; r++)
	g+=is_gap(A->seq_al[r][col]);
      g=(100*g)/A->nseq;
     
      /*enthropy
      for (e=0, r=0; r<26; r++)
	{
	  a0=analyze[r][0];
	  a1=analyze[r][1];
	  t=a0+a1;
	  
	  if (t>0)
	    e+= t/(double)A->nseq*log(t/(double)A->nseq);
	}
      e*=-1;
      */
      e=0;
      if (g>10) continue;
      if (e>10) continue;
      
      if ( strstr ( mode, "SIMTEST"))
	{
	  for (r=0; r<26; r++)
	    {
	      
	      a0=analyze[r][0];
	      a1=analyze[r][1];
	      
	      if (a1)
		{
		  cache[r][col]=a1;
		  cache[26][col]=MAX(cache[26][col],a1);
		}
	    }
	}
      else
	{
	
	  
	  
	  for (r=0; r<26; r++)
	    {
	      
	      a0=analyze[r][0];
	      a1=analyze[r][1];
	      
	      if (!a1 && a0)
		{
		  cache[r][col]=a0;
		  cache[26][col]=MAX(cache[26][col],a0);
		}
	    }
	}
      
      for ( r=0; r<26; r++)analyze[r][0]=analyze[r][1]=0;
      list[col][0]=col;
      list[col][1]=cache[26][col];
    }

  free_int (analyze, -1);

  sort_int_inv (list, 2, 1, 0, A->len_aln-1);
  
  strget_param ( mode, "_MAXN2_", "100000", "%d", &maxn2);

  for ( col=maxn2; col<A->len_aln; col++)
    for ( r=0; r<=26; r++)cache[r][list[col][0]]=0;
 
  free_int (list, -1);
  return cache;
}


int **sar2cache_proba_new ( Alignment *A, int *ns,int **ls, Alignment *S, char *compound, char *mode)
{
  int col, s, seq,ms,mseq, res, mres, res1, n,maxn1, maxn2,maxn3, t, ci, a,w;

  int **list;

  int N1msa,N1sar, N, N11, N10, N01,N00, SCORE, COL_INDEX, RES;
  int nfield=0;
  int value;
  

  int **cache;
  static int **sim;
  int sw_thr;
  float zscore;
  
  RES=nfield++;COL_INDEX=nfield++;N1msa=nfield++;N1sar=nfield++;N=nfield++;N11=nfield++;N10=nfield++;N01=nfield++;N00=nfield++;SCORE=nfield++;
  ci=name_is_in_list ( compound, S->name, S->nseq, -1);
  cache=declare_int (27, A->len_aln);
  
  strget_param ( mode, "_SWTHR_", "30", "%d", &sw_thr);
  strget_param (mode, "_ZSCORE_","0", "%f", &zscore);
  
 
  if (!sim)sim=aln2sim_mat(A, "idmat");
  for ( ms=0; ms<ns[0]; ms++)
    {
      mseq=ls[0][ms];
      if ( S->seq_al[ci][mseq]!='I')continue;

      list=declare_int (A->len_aln+1, nfield);
      for (t=0,n=0, col=0; col< A->len_aln; col++)
	{
	  int same_res;
	  
	  mres=tolower(A->seq_al[mseq][col]);
	  if ( is_gap(mres))continue;
	  
	  list[col][RES]=mres;
	  list[col][COL_INDEX]=col;

	  for ( s=0; s<ns[0]; s++)
	    {
	      seq=ls[0][s];
	      res=tolower(A->seq_al[seq][col]);
	      if (is_gap(res))continue;
	      w=sim[seq][mseq];w=(mres==res)?100-w:w;
	      if (w<sw_thr)w=0;
	      same_res=(res==mres)?1:0;
	      
	      list[col][N]+=w;
	      
	      if (S->seq_al[ci][seq]=='I' && same_res)list[col][N11]+=w;
	      else if (S->seq_al[ci][seq]=='I' && same_res)list[col][N10]+=w;
	      else if (S->seq_al[ci][seq]=='O' && same_res)list[col][N01]+=w;
	      else if (S->seq_al[ci][seq]=='O' && same_res)list[col][N00]+=w;
	      
	      if ( S->seq_al[ci][seq]=='I')list[col][N1sar]+=w;
	      if ( same_res)list[col][N1msa]+=w;
	      
	    }

	  list[col][SCORE]=(int)evaluate_sar_score1 (list[col][N], list[col][N11], list[col][N1msa], list[col][N1sar]);
	  
	}
      strget_param ( mode, "_MAXN1_", "5", "%d", &maxn1);
      sort_int_inv (list,nfield,SCORE,0,A->len_aln-1);
      for (a=0; a<maxn1; a++)
	{
	  col=list[a][COL_INDEX];
	  res1=list[a][RES];
	  value=list[a][SCORE];

	  if ( res1!=0)
	    {
	      cache[res1-'a'][col]+= value;
	      cache[26][col]+=value;
	    }
	}
      free_int (list, -1);
    }

  /*Filter Columns*/
  list=declare_int (A->len_aln+1, nfield);
  for ( col=0; col< A->len_aln; col++)
    {
      list[col][COL_INDEX]=col;
      list[col][SCORE]=cache[26][col];
    }
  /*Keep The N Best Columns*/
  if ( zscore!=0)
    {
      double sum=0, sum2=0, z;
      int n=0;
      for (a=0; a< A->len_aln; a++)
	{
	  if ( list[a][SCORE]>0)
	    {
	      sum+=list[a][SCORE];
	      sum2+=list[a][SCORE]*list[a][SCORE];
	      n++;
	    }
	}
      for (a=0; a<A->len_aln; a++)
	{
	  if ( list[a][SCORE]>0)
	    {
	      z=return_z_score (list[a][SCORE], sum, sum2,n);
	      if ((float)z<zscore)
		{  
		  col=list[a][COL_INDEX];
		  for (s=0; s<27; s++)
		    cache [s][col]=0;
		}
	      else
		{
		  fprintf ( stdout, "\nZSCORE: KEEP COL %d SCORE: %f SCORE: %d\n", list[a][COL_INDEX], (float)z, list[a][SCORE]);
		}
	    }
	}
    }
  else
    {
      sort_int_inv (list,nfield,SCORE,0,A->len_aln-1);
      strget_param ( mode, "_MAXN2_", "100000", "%d", &maxn2);
      
      for (a=maxn2;a<A->len_aln; a++)
	{
	  col=list[a][COL_INDEX];
	  for (s=0; s<27; s++)
	    cache [s][col]=0;
	}
    }

  /*Get Rid of the N best Columns*/;
  strget_param ( mode, "_MAXN3_", "0", "%d", &maxn3);
  
  for (a=0; a<maxn3;a++)
    {
      col=list[a][COL_INDEX];
      for (s=0; s<27; s++)
	cache [s][col]=0;
    }
  return cache;    
}
int **sar2cache_adriana ( Alignment *A, int *ns,int **ls, Alignment *S, char *compound, char *mode)
{
  
  int col,maxn1, s, seq,ms,mseq, res, mres,res1, n, t, ci, a;
  float quant=0;
  int **list;


  int **cache;
  ci=name_is_in_list ( compound, S->name, S->nseq, -1);
  cache=declare_int (A->nseq, A->len_aln);
    
  
  for ( ms=0; ms<ns[0]; ms++)
    {
      mseq=ls[0][ms];
      if ( S->seq_al[ci][mseq]!='I')continue;

      list=declare_int (A->len_aln+1, 5);
      for (t=0,n=0, col=0; col< A->len_aln; col++)
	{
	  mres=tolower(A->seq_al[mseq][col]);
	  list[col][0]=mres;
	  list[col][1]=col;

	  if ( is_gap(mres))continue;
	  for ( s=0; s<ns[0]; s++)
	    {
	      seq=ls[0][s];
	      res=tolower(A->seq_al[seq][col]);
	      if (is_gap(res))continue;
	      
	      if (S->seq_al[ci][seq]=='I' && res==mres)list[col][3]++;
	      if (res==mres)list[col][2]++;
	    }
	}
      
      sort_int_inv (list,5,3,0,A->len_aln-1);
      
      strget_param ( mode, "_MAXN1_", "5", "%d", &maxn1);
      strget_param ( mode, "_QUANT_", "0.95", "%f", &quant);
      
      n=quantile_rank ( list, 3, A->len_aln,quant);
      sort_int (list, 5, 2, 0, n-1);
      
      for (a=0; a<maxn1; a++)
	{
	 
	  col=list[a][1];
	  res1=list[a][0];
	  cache[mseq][col]=list[a][3];
	}
      free_int (list, -1);
     
    }
  return cache;
}
int **sar2cache_proba2 ( Alignment *A, int *ns,int **ls, Alignment *S, char *compound, char *mode)
{
  int col, s, seq,ms,mseq, res, mres,n,maxn1, t, ci, a,b;
  int COL, SCORE;
  
  float quant=0;
  int **list;

  float T1, T2, T3, T4;

  int **cache;
  cache=declare_int ( A->nseq, A->len_aln);
  ci=name_is_in_list ( compound, S->name, S->nseq, -1);
      
  strget_param ( mode, "_FILTER1_", "0"   , "%f", &T1);
  strget_param ( mode, "_FILTER2_", "1000000", "%f", &T2);
  strget_param ( mode, "_FILTER3_", "0"   , "%f", &T3);
  strget_param ( mode, "_FILTER4_", "1000000", "%f", &T4);
  
  list=declare_int (A->len_aln+1,A->nseq+2);
  SCORE=A->nseq;
  COL=A->nseq+1;
  
  for ( ms=0; ms<ns[0]; ms++)
    {
      mseq=ls[0][ms];
      if ( S->seq_al[ci][mseq]!='I')continue;

      for (t=0,n=0, col=0; col< A->len_aln; col++)
	{
	  int N11=0,N10=0,N01=0,N00=0,N1sar=0,N1msa=0,N=0;
	  	  
	  mres=tolower(A->seq_al[mseq][col]);
	  if ( is_gap(mres))continue;
	  for ( s=0; s<ns[0]; s++)
	    {
	      seq=ls[0][s];
	      res=tolower(A->seq_al[seq][col]);
	      if (is_gap(res))continue;
	      
	      N++;
	      if (S->seq_al[ci][seq]=='I' && res==mres)N11++;
	      else if (S->seq_al[ci][seq]=='I' && res!=mres)N10++;
	      else if (S->seq_al[ci][seq]=='O' && res==mres)N01++;
	      else if (S->seq_al[ci][seq]=='O' && res!=mres)N00++;

	      if ( S->seq_al[ci][seq]=='I')N1sar++;
	      if ( res==mres)N1msa++;
	    }
	  list[col][mseq]=(int)evaluate_sar_score1 (N,N11,N1msa,N1sar);
	  list[col][SCORE]+=list[col][mseq];
	  list[col][COL]=col;
	}
    }

  strget_param ( mode, "_MAXN1_", "5", "%d", &maxn1);
  strget_param ( mode, "_QUANT_", "0.95", "%f", &quant);
  sort_int_inv (list,A->nseq+2,SCORE, 0, A->len_aln-1);
  n=quantile_rank ( list,A->nseq, A->len_aln,quant);
  n=5;
  

  for (a=0; a<n; a++)
    {
      int value;
      
      col=list[a][COL];
      for ( b=0; b<A->nseq; b++)
	{
	  value=list[col][b];
	  if ( value>T1 && value<T2){cache[b][col]= value;}  
	}
    }
  
  free_int (list, -1);
  return cache;
}

  
	  

/************************************************************************************/
/*                ALIGNMENT ANALYZE     : SAR                                            */
/************************************************************************************/
int aln2jack_group3 (Alignment *A,char *comp, int **l1, int *nl1, int **l2, int *nl2)
{
  int **seq_list, **sar_list, nsar=0, nseq=0;
  int a, b, mid;

  vsrand (0);
  sar_list=declare_int (A->nseq, 2);
  seq_list=declare_int (A->nseq, 2);
  for (a=0; a< A->nseq; a++)
    {
      if (comp[a]=='I')
	{
	  sar_list[nsar][0]=a;
	  sar_list[nsar][1]=rand()%100000;
	  nsar++;
	}
      else
	{
	  seq_list[nseq][0]=a;
	  seq_list[nseq][1]=rand()%100000;
	  nseq++;
	}
    }
  
  
  l1[0]=vcalloc (A->nseq, sizeof (int));
  l2[0]=vcalloc (A->nseq, sizeof (int));
  nl1[0]=nl2[0]=0;
  
  sort_int (seq_list, 2, 1, 0,nseq-1);
  sort_int (sar_list, 2, 1, 0,nsar-1);
  mid=nsar/2;
  for (a=0; a<mid; a++)
    {
      l1[0][nl1[0]++]=sar_list[a][0];
    }
  for (a=0,b=mid; b<nsar; b++, a++)
    {
      l2[0][nl2[0]++]=sar_list[b][0];
    }

  mid=nseq/2;
  for (a=0; a<mid; a++)
    {
      l1[0][nl1[0]++]=seq_list[a][0];
    }
  for (a=0,b=mid; b<nseq; b++, a++)
    {
      l2[0][nl2[0]++]=seq_list[b][0];
    }

  
  free_int (seq_list, -1);
  free_int (sar_list, -1);
  return 1;
}

int aln2jack_group2 (Alignment *A, int seq, int **l1, int *nl1, int **l2, int *nl2)
{
  int **list;
  int a, b, mid;
  

  list=declare_int (A->nseq, 2);
  l1[0]=vcalloc (A->nseq, sizeof (int));
  l2[0]=vcalloc (A->nseq, sizeof (int));
  nl1[0]=nl2[0];
  
  vsrand (0);
  for ( a=0; a< A->nseq; a++)
    {
      list[a][0]=a;
      list[a][1]=rand()%100000;
    }
  sort_int (list, 2, 1, 0,A->nseq-1);
  mid=A->nseq/2;
  for (a=0; a<mid; a++)
    {
      l1[0][nl1[0]++]=list[a][0];
    }
  for (a=0,b=mid; b<A->nseq; b++, a++)
    {
      l2[0][nl2[0]++]=list[b][0];
    }

  free_int (list, -1);
  return 1;
}
int aln2jack_group1 (Alignment *A, int seq, int **l1, int *nl1, int **l2, int *nl2)
{
  int **sim;
  int **list;
  int a, mid;
  
  list=declare_int ( A->nseq, 3);
  l1[0]=vcalloc (A->nseq, sizeof (int));
  l2[0]=vcalloc (A->nseq, sizeof (int));
  nl1[0]=nl2[0];
  
  sim=aln2sim_mat (A, "idmat");
  for ( a=0; a< A->nseq; a++)
    {
      list[a][0]=seq;
      list[a][1]=a;
      list[a][2]=(a==seq)?100:sim[seq][a];
    }
  sort_int_inv (list, 3, 2, 0, A->nseq-1);
  fprintf ( stderr, "\nJacknife fromsequence %s [%d]\n", A->name[seq], seq);
  mid=A->nseq/2;
  for (a=0; a< mid; a++)
    l1[0][nl1[0]++]=list[a][1];
  for (a=mid; a<A->nseq; a++)
    l2[0][nl2[0]++]=list[a][1];
  return 1;
}
  
      
int sarset2subsarset ( Alignment *A, Alignment *S, Alignment **subA, Alignment **subS, Alignment *SUB)
{
  Alignment *rotS, *intS;
  int a,b, *list, nl;
  
  list=vcalloc ( SUB->nseq, sizeof (int));
  for (nl=0,a=0; a<SUB->nseq; a++)
    {
      b=name_is_in_list(SUB->name[a], A->name, A->nseq, 100);
      if ( b!=-1)list[nl++]=b;
    }

  subA[0]=extract_sub_aln (A, nl, list);
  rotS=rotate_aln (S, NULL);
  intS=extract_sub_aln (rotS, nl, list);
    
  subS[0]=rotate_aln (intS, NULL);

  for ( a=0; a<S->nseq; a++) sprintf ( (subS[0])->name[a], "%s", S->name[a]);
  
  
  return 0;
}

int ***simple_sar_analyze_vot ( Alignment *A, Alignment *SAR, char *mode)
{
  int a, b, c, d;
  int res1, res2, sar1, sar2;
  float s;
  int **sim;
  static float ***result;
  static int ***iresult;
  if (!result)
    {
    result=declare_arrayN (3,sizeof (float),SAR->nseq, A->len_aln,3);
    iresult=declare_arrayN (3,sizeof (int),SAR->nseq, A->len_aln,3);
    }

  sim=aln2sim_mat (A, "idmat");
  
  
  for (a=0; a<SAR->nseq; a++)
    for (b=0; b<A->len_aln; b++)
      result[a][b][0]=1;
 
  for ( a=0; a< SAR->nseq; a++)
    for ( b=0; b<A->nseq-1; b++)
      for ( c=b+1; c< A->nseq; c++)
	for ( d=0; d<A->len_aln; d++)
	  {
	    res1=A->seq_al[b][d];
	    res2=A->seq_al[c][d];

	    sar1=(SAR->seq_al[a][b]=='I')?1:0;
	    sar2=(SAR->seq_al[a][c]=='I')?1:0;
	    
	    s=sim[b][c];
	    
	    
	    
	    
	    if ( sar1!=sar2 && res1!=res2)
	      result[a][d][0]*=(1/(100-s));
	    
	    else if ( sar1==sar2 && sar1==1 && res1==res2)
	      result[a][d][0]*=1/s;
	    
	    
	    

	    /*
	    else if ( sar1==sar2 && res1==res2)result[a][d][0]+=(100-s)*(100-s);
	    else if ( sar1==sar2 && res1!=res2)result[a][d][0]-=s*s;
	    else if ( sar1!=sar2 && res1==res2)result[a][d][0]-=(100-s)*(100-s);
	    */
	    
	    result[a][d][1]='a';
	  }
  for ( a=0; a<SAR->nseq; a++)
    for ( b=0; b<A->len_aln; b++)
      {
	fprintf ( stderr, "\n%f", result[a][b][0]);
	iresult[a][b][0]=100*log(1-result[a][b][0]);
      }
  return iresult;
}
int display_simple_sar_analyze_pair_col (Alignment *A, Alignment *SAR, char *mode)
{
  int **r;
  int a, b, n, do_tm;

 
  Alignment *rA;
  int *nI;

  
  strget_param (mode, "_TM_", "0", "%d", &do_tm);
  r=simple_sar_analyze_pair_col (A, SAR,mode);
  rA=rotate_aln (A, NULL);
  n=0;
  
  nI=vcalloc ( SAR->nseq, sizeof (int));
  for (a=0; a< SAR->nseq; a++)
    for (b=0; b<SAR->len_aln; b++) nI[a]+=(SAR->seq_al[a][b]=='I')?1:0;
  
    
  
  while ( r[n][0]!=-1)
    {
      if (r[n][3]>0)
	{fprintf ( stdout,  "COMP S: %3d %3d %s %20s %2d #\n", r[n][3],0,SAR->seq_al[r[n][0]], SAR->name[r[n][0]], nI[r[n][0]]);
	  fprintf ( stdout, "SEQ1 S: %3d %3d %s %20s %2d #\n", r[n][3],r[n][1],rA->seq_al[r[n][1]], SAR->name[r[n][0]],nI[r[n][0]]);
	  fprintf ( stdout, "SEQ2 S: %3d %3d %s %20s %2d #\n\n", r[n][3],r[n][2],rA->seq_al[r[n][2]], SAR->name[r[n][0]],nI[r[n][0]]);
	}
      n++;
    }
  return 0;
}
int display_simple_sar_analyze_col (Alignment *A, Alignment *SAR, char *mode)
{
  int ***result, **r2, **r3, **r4, **aa;
  int a, b, c, n;
  char *cons;
  int threshold=20;
  int do_tm;
  strget_param (mode, "_TM_", "0", "%d", &do_tm);
  result=simple_sar_analyze_col (A, SAR,mode);
  r2=declare_int (A->len_aln*SAR->nseq, 5);
  r3=declare_int (A->len_aln+1, 5);
  r4=declare_int (A->len_aln+1, SAR->nseq+1);
  aa=declare_int (2, 256);
  cons=vcalloc (A->len_aln+1, sizeof (char));
  for (a=0; a<A->len_aln; a++){r3[a][0]=a;cons[a]='A';}
  
  
  
  for (n=0,a=0; a< SAR->nseq; a++)
    {
      double sum, sum2;
      for (sum=0, sum2=0,b=0; b<A->len_aln; b++)
	{
	  sum+=result[a][b][0];
	  sum2+=result[a][b][0]*result[a][b][0];
	}
      
      for (b=0; b<A->len_aln; b++, n++)
	{
	  r2[n][0]=a;//compound
	  r2[n][1]=b;//pos
	  r2[n][2]=result[a][b][1]; //AA
	  r2[n][3]=result[a][b][0]; //Score
	  r2[n][4]=result[a][b][2];  //(int)10*return_z_score ((double)result[a][b][0], sum, sum2, A->len_aln); //ZScore
      }
    }
  sort_int (r2,5, 3, 0, n-1);//sort on Score (3rd field)
  for ( a=0; a< n; a++)
    {
      int comp, pos, bad;

      comp=r2[a][0];
      pos=r2[a][1];
      fprintf ( stdout, "SEQ %5d %5d %5d %s ",r2[a][1]+1,r2[a][3], r2[a][4], (do_tm)?alnpos2hmmtop_pred (A, NULL, r2[a][1], SHORT):"x");
      for (c=0; c<A->nseq; c++)fprintf (stdout, "%c", A->seq_al[c][r2[a][1]]);
     
      
      
      bad=0;
      for (c=0; c< A->nseq; c++)
	{
	  int activity, res;
	  
	  activity=SAR->seq_al[comp][c];
	  res=A->seq_al[c][pos];
	  
	  if (activity=='O')aa[0][res]++;
	  if (activity=='I')aa[1][res]++;
	}
            
      for (c=0; c< A->nseq; c++)
	{
	  int activity, res;
	  activity=SAR->seq_al[comp][c];
	  res=A->seq_al[c][pos];
	  bad+=(aa[0][res] && aa[1][res])?1:0;
	  aa[0][res]=aa[1][res]=0;
	} 
      fprintf ( stdout, " %20s %d |\nCOM %5d %5d %5d %s %s %20s %d |\n\n", SAR->name[r2[a][0]],bad,r2[a][1]+1,r2[a][3],r2[a][4], (do_tm)?alnpos2hmmtop_pred (A, NULL, r2[a][1], SHORT):"x",SAR->seq_al[r2[a][0]], SAR->name[r2[a][0]], bad);

      
      if (r2[a][4]>threshold)
	{
	  cons[r2[a][1]]++;
	  r3[r2[a][1]][1]++;
	  r3[r2[a][1]][2]+=r2[a][3];
	  r4[r2[a][1]][r2[a][0]]=1;
	}
    }
  sort_int (r3, 3,1,0, A->len_aln-1);

  for (a=0; a<A->len_aln; a++)
    {
      if (r3[a][1]>0)
	{
	  fprintf ( stdout, "\nPOS %4d %4d %4d %c ", r3[a][0]+1, r3[a][1], r3[a][2], cons[r3[a][0]]);
	  for (b=0; b<SAR->nseq; b++)fprintf ( stdout, "%d", r4[r3[a][0]][b]);
	  if (do_tm)fprintf ( stdout, " %s",alnpos2hmmtop_pred (A, NULL, r3[a][0], VERBOSE));
      }
    }
  for (a=0; a< A->nseq; a++)fprintf ( stdout, "\n#MSA >%s\n#MSA %s",A->name[a], A->seq_al[a]);
  fprintf ( stdout, "\n#MSA >cons\n#MSA %s", cons);
  
  return 0;
}
int ***    simple_sar_predict     (Alignment *A, Alignment *SAR, char *mode)
{
  //This function estimates the z score of every poition with every compound
  //The best Z-score position is then used for the prediction

  int a, b, c, nts, pos, Rscore,Zscore;
  int ***r;
  int ***pred;
  int **aa;
  

  aa=declare_int (2,256);
  pred=declare_arrayN (3, sizeof (int),SAR->nseq, A->nseq, 5);
  

  r=simple_sar_analyze_col (A, SAR, mode);
  nts=SAR->len_aln; //number_trainning_sequences;
  
  for (a=0; a<SAR->nseq; a++)
    {
      sort_int (r[a],4, 2, 0, A->len_aln-1);
   
      pos=r[a][A->len_aln-1][3]; //Best Position
      Zscore=r[a][A->len_aln-1][2]; //Best Z-Score
      Rscore=r[a][A->len_aln-1][0]; //Best Z-Score
     
     
      for (c=0; c<nts; c++)
	{

	  if (SAR->seq_al[a][c]=='I')aa[1][(int)A->seq_al[c][pos]]++;//Build Positive Alphabet for Compound a
	  if (SAR->seq_al[a][c]=='O')aa[0][(int)A->seq_al[c][pos]]++;//Build Positive Alphabet for Compound a
	}
      for (c=nts; c<A->nseq; c++)
	{
	  pred[a][c][0]=pos;
	  pred[a][c][1]=Zscore;
	  pred[a][c][2]=Rscore;
	  if (aa[1][(int)A->seq_al[c][pos]]>0)
	    {

	      pred[a][c][3]=aa[1][(int)A->seq_al[c][pos]];
	      pred[a][c][4]=aa[0][(int)A->seq_al[c][pos]];
	    }
	}
      for (c=0; c<nts; c++)aa[0][(int)A->seq_al[c][pos]]=aa[1][(int)A->seq_al[c][pos]]=0;
    }

  for ( a=nts; a< A->nseq; a++)
    {
      for ( b=0; b<SAR->nseq; b++)
	{
	  fprintf ( stdout, ">%-25s %-25s Pos %3d ZScore %3d Rscore %3d Activity +: %d -: %d ", A->name [a], SAR->name[b], pred[b][a][0],pred[b][a][1], pred[b][a][2], pred[b][a][3], pred[b][a][4]);
	  if (pred[b][a][4]==0)for (c=0; c<pred[b][a][3]; c++)fprintf ( stdout, "*");
	  fprintf ( stdout, "\n" );
	  for (c=0; c<A->nseq; c++)fprintf ( stdout, "%c", A->seq_al[c][ pred[b][a][0]]);
	  fprintf ( stdout, " %s\n", SAR->name[b]);
	  for (c=0; c<A->nseq-1; c++)fprintf ( stdout, "%c", SAR->seq_al[b][c]);
	  fprintf ( stdout, "  %s\n", SAR->name[b]);	  fprintf ( stdout, "\n");
	  
	}
    }
  return pred;
}
int *pair_seq2seq (int *iseq, char *seq1, char *seq2);
int **simple_sar_analyze_pair_col ( Alignment *inA, Alignment *SAR, char *mode)
{

  int a, b, c, n, n2;
  int *iseq=NULL;
  static int **result, **fresult;
  int sar_mode=1;
  int maxgapratio=0;
  int nresults=10;
  double sum, sum2, score;
  Alignment *A;
  char aa;

  if (!result)
    {
      result=declare_int (inA->len_aln*inA->len_aln,5);
      
      fresult=declare_int (inA->len_aln*nresults*SAR->nseq, 5);
      
    }

  A=rotate_aln (inA, NULL);
  

  for (n2=0,a=0; a<SAR->nseq; a++)
    {

      for (n=0, sum=0, sum2=0,b=0; b<A->nseq-1; b++)
	{
	  for ( c=b+1; c<A->nseq; c++, n++)
	    {

	      iseq=pair_seq2seq (iseq,A->seq_al[b], A->seq_al[c]);
	      if ( sar_mode==1)
		score=sar_vs_iseq1(SAR->seq_al[a],iseq,maxgapratio,NULL,&aa);
	      else if (sar_mode==4)
		score=sar_vs_iseq4(SAR->seq_al[a],iseq,maxgapratio,NULL,&aa);
	      //HERE ("%d", (int)score);
	      result[n][0]=a;//compound;
	      result[n][1]=b; //pos1
	      result[n][2]=c; //pos2
	      result[n][3]=(int)score;
	      
	      sum+=score;
	      sum2+=score*score;
	    }
	}
      for (b=0; b<n; b++)
	result[c][4]=(int)10*return_z_score ((double)result[c][3], sum, sum2,n); //ZScore
      sort_int (result,5,3,0,n-1);
      for ( b=n-nresults; b<n; b++, n2++)
	{

	  for ( c=0; c<5; c++)fresult[n2][c]=result[b][c];
	}
    }
  sort_int (fresult,5,3,0,n2-1);
  fresult[n2][0]=-1; //Terminator;
  
  return fresult;
}

int *pair_seq2seq (int *iseq, char *seq1, char *seq2)
{
  static int **lu;
  int a;
  if (!lu)
    {
      int a, b, n=0;
      lu=declare_int (256,256);
      for ( n=a=0; a<256; a++)
	for (b=0; b<256; b++)
	  lu[a][b]=n++;
    }
  
  if (!iseq)iseq=vcalloc ( strlen (seq1)+1, sizeof (int));
  for ( a=0; a< strlen (seq1); a++)
    {
      if (is_gap(seq1[a]) || is_gap(seq2[a]))iseq[a]=-1;
      else iseq[a]=lu[(int)seq1[a]][(int)seq2[a]];
    }
  
  return iseq;
}
int ***simple_sar_analyze_col ( Alignment *inA, Alignment *SAR, char *mode)
{
  Alignment *A;
  double score=0, best_score=0;
  int best_pos=0;
  int a, b;
 
  static int ***result;
  int **sim;
  char aa;
  int sar_mode=3;
  int maxgapratio=0;
  double sum, sum2;
  

  strget_param (mode, "_METHOD_", "3", "%d", &sar_mode);
  if (!result)
    result=declare_arrayN (3,sizeof (int),SAR->nseq, inA->len_aln,4);
  

  sim=aln2sim_mat (inA, "idmat");
  A=rotate_aln (inA, NULL);
  
  
  for ( a=0; a<SAR->nseq; a++)
    {
      best_pos=best_score=0;
      for ( sum=0, sum2=0,b=0; b<A->nseq; b++)
	{
	
	  if ( sar_mode==1)
	    score=sar_vs_seq1(SAR->seq_al[a], A->seq_al[b],maxgapratio, sim, &aa);
	  else if ( sar_mode==2)
	    score=sar_vs_seq2(SAR->seq_al[a], A->seq_al[b],maxgapratio, sim, &aa);
	  else if (sar_mode ==3)
	    score=sar_vs_seq3(SAR->seq_al[a], A->seq_al[b],maxgapratio, sim, &aa);
	  else if (sar_mode ==4)
	    score=sar_vs_seq4(SAR->seq_al[a], A->seq_al[b],maxgapratio, sim, &aa);
	  else if (sar_mode ==5)
	    score=sar_vs_seq5(SAR->seq_al[a], A->seq_al[b],maxgapratio, sim, &aa);
	  
	  
	  result[a][b][0]+=score*10;
	  result[a][b][1]=aa;
	  result[a][b][3]=b;
	  sum+=result[a][b][0];
	  sum2+=result[a][b][0]*result[a][b][0];
	  
	}
      for ( b=0; b< A->nseq; b++)result[a][b][2]=10*return_z_score ((double)result[a][b][0], sum, sum2, A->nseq); //Score
    }
  
  return result;

 }
int *seq2iseq ( char *seq);
double sar_vs_seq4 ( char *sar, char *seq, float gl, int **sim, char *best_aa)
{
  
  return sar_vs_iseq4 (sar, seq2iseq(seq), gl, sim, best_aa);
}
double sar_vs_seq1 ( char *sar, char *seq, float gl, int **sim, char *best_aa)
{
  
  return sar_vs_iseq1 (sar, seq2iseq(seq), gl, sim, best_aa);
}

int *seq2iseq ( char *seq)
{
 static int *iseq, clen;
  int a;
  
  if (!iseq || clen<strlen (seq))
    {
      clen=strlen (seq);
      vfree (iseq);
      iseq=vcalloc (clen+1, sizeof (int));
    }
  for (a=0; a<strlen (seq); a++)
    {
      iseq[a]=seq[a];
    }
  return iseq;
}

double sar_vs_iseq1 ( char *sar, int *seq, float gl, int **sim, char *best_aa)
{
  double score=0, return_score=0;
  int RN,N11, Nmsa, Nsar, N, N10, N01, N00;
  int a, b, r, s, res, res1;
  double Ng=0;
  static int *aa, *aal, naa;
  
  /*measure the E-Value for every amino acid. Returns the best one*/


  
  N=strlen (sar);
  for (a=0; a<N; a++)
    Ng+=(is_gap(seq[a])|| seq[a]==-1);
  Ng/=N;
  
  if (Ng>gl) return 0;
  
  if (!aa)
    {
      aa=vcalloc (256*256, sizeof(int));
      aal=vcalloc (N, sizeof (int));
    }
  naa=0;
  for ( a=0; a<N; a++)
    {
      if (!aa[seq[a]])
	{
	  aal[naa++]=seq[a];
	  aa[seq[a]]=1;
	}
    }
  
  best_aa[0]='-';
  for (a=0; a<naa; a++)
    {

      res=aal[a];
  
      if (res==-1 || !aa[res]);
      else
	{
	  RN=Nmsa=Nsar=N11=N10=N01=N00=0;
	  
	  for (b=0; b<N; b++)
	    {
	      	      
	      res1=seq[b];
	      if (res1=='-' || res1==-1)r=0;
	      else r=(res1==res)?1:0;
	      
	     
	      s=(sar[b]=='I')?1:0;
	      
	      Nmsa+=r; Nsar+=s;
	      N11+=(r && s)?1:0;
	      N01+=(!r &&s)?1:0;
	      N10+=(r && !s)?1:0;
	      N00+=(!r && !s)?1:0;
	      RN++;
	    }
	  
	  if (N11)
	    {
	      score=evaluate_sar_score1 ( RN, N11, Nmsa, Nsar);
	      
	    }
	  else
	    {
	      score=0;
	    }
	  
	  if ( score>return_score)
	    {
	      best_aa[0]=res;
	      return_score=score;
	    }
	}
    }
  for ( a=0; a<N; a++)aa[seq[a]]=0;
  
  return return_score;
}
double sar_vs_seq5 ( char *sar, char *seq, float gl, int **sim, char *best_aa)
{
  double score=0;
  int N11, Nmsa, Nsar, N, N10, N01, N00;
  int a, b, r, s;
  double Ng=0;
  int *aa;

  //measure the E-Value if all the 1AA are considered like alphabet 1*/
  //This function is testing ONE Compound (sar) against one column of MSA (seq) containing N sequences
  //Returns 0 if the active and inactive alphabet overlap  
  N=strlen (sar);
  //N is the nmber of sequences

  for (a=0; a<N; a++)
    Ng+=is_gap(seq[a]);
  Ng/=N;
  
  if (Ng>gl) return 0;

  //Identify all the AA associated with a I (Positive alphabet)
  aa=vcalloc ( 256, sizeof (int));  
  for (b=0; b<N; b++)
    {
      s=(sar[b]=='I')?1:0;
      if (s)aa[(int)seq[b]]=1;
    }
  N11=N10=N01=N00=Nmsa=Nsar=0;


  for (b=0; b<N; b++)
    {
      
      r=aa[(int)seq[b]];
      s=(sar[b]=='I')?1:0;
	  
      Nmsa+=r; Nsar+=s;
      N11+=(r && s)?1:0;
      N01+=(!r &&s)?1:0;
      N10+=(r && !s)?1:0;
      N00+=(!r && !s)?1:0;
    }
  if (N10>=1 || N01>=1) return 0;
  if (N11)
    {
      score=evaluate_sar_score1 ( N, N11, Nmsa, Nsar);
    }
  else score=0;
  
  vfree (aa);
  return score;
  
}


  
double sar_vs_iseq4 ( char *sar, int *seq, float gl, int **sim, char *best_aa)
{
  int N, Ni, No;
  int a, b,c, r, s;
  double Ng=0;
  static int **aa;
  
  /*Correlation between AA conservation and Activity*/
  
  N=strlen (sar);
  for (a=0; a<N; a++)
    Ng+=(is_gap(seq[a]) || seq[a]==-1)?1:0;
  Ng/=N;
  if (gl<1)Ng*=100;
  
  if (Ng>gl) return 0;
  

  if (!aa)aa=declare_int(2,257*257);
  for (No=Ni=b=0; b<N; b++)
    {

      s=(sar[b]=='I')?1:0;
      if (s){aa[1][seq[b]]++;Ni++;}
      else  {aa[0][seq[b]]++;No++;}
     }
  for (r=0,b=0; b<N; b++)
    {
      if (aa[1][(int) seq[b]]==Ni)
	{
	  r=(No-aa[0][(int)seq [b]])*100/No;
	  break;
	}
      
    }
  for (c=0; c<N; c++)aa[0][seq[c]]=aa[1][seq[c]]=0;
  
  return r;
}

double sar_vs_seq3 ( char *sar, char *seq, float gl, int **sim, char *best_aa)
{
  double score=0;
  int N11, Nmsa, Nsar, N, N10, N01, N00;
  int a, b, r, s;
  double Ng=0;
  int *aa;

  //measure the E-Value if all the 1AA are considered like alphabet 1*/
  //This function is testing ONE Compound (sar) against one column of MSA (seq) containing N sequences
    
  N=strlen (sar);
  //N is the nmber of sequences

  for (a=0; a<N; a++)
    Ng+=is_gap(seq[a]);
  Ng/=N;
  
  if (Ng>gl) return 0;

  //Identify all the AA associated with a I (Positive alphabet)
  aa=vcalloc ( 256, sizeof (int));  
  for (b=0; b<N; b++)
    {
      s=(sar[b]=='I')?1:0;
      if (s)aa[(int)seq[b]]=1;
    }
  N11=N10=N01=N00=Nmsa=Nsar=0;


  for (b=0; b<N; b++)
    {
      
      r=aa[(int)seq[b]];
      s=(sar[b]=='I')?1:0;
	  
      Nmsa+=r; Nsar+=s;
      N11+=(r && s)?1:0;
      N01+=(!r &&s)?1:0;
      N10+=(r && !s)?1:0;
      N00+=(!r && !s)?1:0;
    }
  
  if (N11)
    {
      score=evaluate_sar_score1 ( N, N11, Nmsa, Nsar);
    }
  else score=0;
  
  vfree (aa);
  return score;
  
}

double sar_vs_seq2 ( char *sar, char *seq, float gl, int **sim_mat, char *best_aa)
{
  double score=0, return_score=0;
  int L,N11, Nmsa, Nsar,N10, N01, N;
  int a, b,c,d, r1, s1,r2, s2, res;
  double Ng=0;
  int sim, diff, w;
  char string[5];

  /*Weighted E-Value Similarity*/
  L=strlen (sar);
  for (a=0; a<L; a++)
    Ng+=is_gap(seq[a]);
  Ng/=L;
  
  if (Ng>gl) return 0;
  for (a=0; a<26; a++)
    {

      N=Nmsa=Nsar=N11=N10=N01=0;
      res='a'+a;
      for (d=0,b=0; b<L; b++)d+=((tolower(seq[b]))==res)?1:0;
      if ( d==0) continue;
      
      for (b=0; b<L; b++)
	{
	  r1=(tolower(seq[b])==res)?1:0;
	  s1=(sar[b]=='I')?1:0;
	  for ( c=0; c<L; c++)
	    {
	      r2=(tolower(seq[c])==res)?1:0;
	      s2=(sar[c]=='I')?1:0;
	    
	      sprintf ( string, "%d%d%d%d", r1,s1, r2, s2);
	      sim= sim_mat[b][c]/10;
	      diff=10-sim;
	      
	      if (strm (string, "0000"))      {w=diff;N+=2*w;}
	      else if ( strm (string, "0011")){w=sim ;N+=2*w ; N11+=w  ;N10+=0   ;N01+=w   ;Nmsa+=w   ;Nsar+=w;}
	      else if ( strm (string, "1010")){w=diff;N+=2*w ; N11+=0  ;N10+=2*w ;N01+=0   ;Nmsa+=2*w ;Nsar+=0;}
	      else if ( strm (string, "0101")){w=diff;N+=2*w;  N11+=0  ;N10+=0   ;N01+=2*w ;Nmsa+=0   ;Nsar+=2*w;}
	      else if ( strm (string, "1111")){w=diff;N+=2*w;  N11+=2*w;N10+=0   ;N01+=0   ;Nmsa+=2*w ;Nsar+=2*w;}
	      else if ( strm (string, "1001")){w=sim; N+=2*w;  N11+=0  ;N10+=w   ;N01+=w   ;Nmsa+=w;Nsar+=w;}
	      else if ( strm (string, "0110")){w=sim; N+=2*w;  N11+=0  ;N10+=w   ;N01+=w   ;Nmsa+=w;Nsar+=w;}
	    }
	}
      if (N11)
	{
	 
	  score=evaluate_sar_score1 ( N, N11, Nmsa, Nsar);
	}
      return_score=MAX(return_score, score);
    }
  if ( return_score <0)fprintf ( stderr, "\n%.2f", return_score);
  return return_score;
}
  
float get_sar_sim (char *seq1, char *seq2)
{
  int a, l, s, r;
  int n11=0, n10=0, n01=0, n00=0;
  

  l=strlen (seq1);
  for ( a=0; a<l; a++)
    {
      s=(seq1[a]=='O')?0:1;
      r=(seq2[a]=='O')?0:1;

      n00+=(!s && !r)?1:0;
      n11+=(s && r)?1:0;
      n01+=(!s && r)?1:0;
      n10+=(s && !r)?1:0;
    }
  if ( n11==0) return 0;
  else return ((float)(n11)*100)/(float)(n11+n10+n01);
}
float get_sar_sim2 (char *seq1, char *seq2)
{
  int a, l, s, r;
  int n11=0, n10=0, n01=0, n00=0, Ns=0, Nr=0;
  float score1;

  l=strlen (seq1);
  for ( a=0; a<l; a++)
    {
      s=(seq1[a]=='O')?0:1;
      r=(seq2[a]=='O')?0:1;
      
      Ns+=s; 
      Nr+=r;
      
      n00+=(!s && !r)?1:0;
      n11+=(s && r)?1:0;
      n01+=(!s && r)?1:0;
      n10+=(s && !r)?1:0;
    }
  if ( n11==0) return 0;
  



  score1=evaluate_sar_score1 (l, n11, Ns, Nr);
  
  return score1;
}	  
float sar_aln2cor (Alignment *A);
int   sar_aln2ev  (Alignment *A);

int sarseq2ev ( char *s1, char *s2, int MODE);
char*  sarseq2anti_sarseq (char *seq_in, char *seq_out);
Sequence * compare_sar_sequence( Sequence *S1, Sequence *S2, int depth)
{
  int a, b, c, d, e, f;
  char *s1, *s2, *s3,*s4,*s5, *n1, *n2, *n3, *n4, *n5;
  int max=0;;
  int max_depth=5;
  int score;
  char *a1, *a0, *a2, *a3, *a4;
  Alignment *A;

  if ( depth>max_depth)
    {
      printf_exit (EXIT_FAILURE, stderr,"maximum depth: %d", max_depth);
    }
  if ( depth==0) depth=2;
  A=declare_aln2 (strlen (S1->seq[0]),depth);
  a0=A->seq_al[0];
  a1=A->seq_al[1];
  A->len_aln=strlen (S1->seq[0]);
  
  for (a=0; a< S1->nseq; a++)
    for ( b=0; b<S2->nseq; b++)
      {
	A->nseq=2;
	sprintf (a0, "%s", S1->seq[a]);
	sprintf (a1, "%s", S2->seq[b]);
	
	if ( strlen (a0)!=strlen (a1))
	  {
	    add_warning (stderr, "WARNING %s (%d) and %s (%d) do not have the same length", S1->name[a], strlen (S1->seq[a]), S2->name[b], strlen (S2->seq[b]));
	    exit (0);
	  }
	
	fprintf ( stdout, ">2 %15s %15s CORR: %.3f EVAL: %5d\n",S1->name[a], S2->name[b],  sar_aln2cor (A), sar_aln2ev (A));
	sarseq2anti_sarseq (S1->seq[a],a0);
	fprintf ( stdout, ">2 %15s %15s ANTI: %.3f EVAL: %5d\n", S1->name[a], S2->name[b], sar_aln2cor (A), sar_aln2ev (A));
	if ( depth >=3)
	  {
	    A->nseq=3;
	    a2=A->seq_al[2];
	    for (c=b+1; c<S2->nseq; c++)
	      {
		sprintf (a0, "%s", S1->seq[a]);
		sprintf (a1, "%s", S2->seq[b]);
		sprintf (a2, "%s", S2->seq[c]);
		fprintf ( stdout, ">2 %15s %15s %15s CORR: %.3f EVAL: %5d\n",S1->name[a], S2->name[b],S2->name[c],  sar_aln2cor (A), sar_aln2ev (A));
		sarseq2anti_sarseq (S1->seq[a],a0);
		fprintf ( stdout, ">2 %15s %15s ANTI: %.3f EVAL: %5d\n", S1->name[a], S2->name[b], S2->name[c],sar_aln2cor (A), sar_aln2ev (A));
	      }
	    if ( depth>=4)
	      {
		A->nseq=4;
		a3=A->seq_al[2];
		for (d=c+1; d<S2->nseq; d++)
		  {
		    sprintf (a0, "%s", S1->seq[a]);
		    sprintf (a1, "%s", S2->seq[b]);
		    sprintf (a2, "%s", S2->seq[c]);
		    sprintf (a3, "%s", S2->seq[d]);
		    
		    fprintf ( stdout, ">2 %15s %15s %15s %15s CORR: %.3f EVAL: %5d\n",S1->name[a], S2->name[b],S2->name[c],S2->name[d],  sar_aln2cor (A), sar_aln2ev (A));
		    sarseq2anti_sarseq (S1->seq[a],a0);
		    fprintf ( stdout, ">2 %15s %15s %15s %15s ANTI: %.3f EVAL: %5d\n", S1->name[a], S2->name[b], S2->name[c],S2->name[d],sar_aln2cor (A), sar_aln2ev (A));
		  }
		if (depth>=5)
		  {
		    A->nseq=5;
		    a4=A->seq_al[3];
		    for (e=d+1; e<S2->nseq; e++)
		      {
			sprintf (a0, "%s", S1->seq[a]);
			sprintf (a1, "%s", S2->seq[b]);
			sprintf (a2, "%s", S2->seq[c]);
			sprintf (a3, "%s", S2->seq[d]);
			sprintf (a4, "%s", S2->seq[d]);
		      }
		  }
	      }
	  }
      }

  return S1;
}
char*  sarseq2anti_sarseq (char *seq_in, char *seq_out)
{
  int a;
  if (!seq_out)seq_out=vcalloc (strlen (seq_in)+1, sizeof (char));
  for (a=0; a<strlen (seq_in); a++)seq_out[a]=(seq_in[a]=='I')?'O':'I';
  return seq_out;
}
float sar_aln2cor (Alignment *A)
{
  float n1, n11, tot_cor;
  int a, b, c,n;
  
  tot_cor=0;
  for (n=0,a=0; a<A->nseq-1; a++)
    for (b=a+1; b<A->nseq; b++)
      {
	for (n11=n1=0,c=0; c<A->len_aln; c++)
	  {
	    n11+=(A->seq_al[a][c]=='I' && A->seq_al[b][c]=='I');
	    n1+= (A->seq_al[a][c]=='I' || A->seq_al[b][c]=='I');
	  }
	tot_cor+=(n1==0)?0:n11/n1;
	n++;
      }
  tot_cor/=n;
  return tot_cor;
}
int sarseq_pair2ev  ( char *s1, char *s2,int mode);
int sar_aln2ev (Alignment *A)
{
  float n1, n11;
  int a, b, c, tot=0, n=0;
  
  tot=0;
  for (a=0; a<A->nseq-1; a++)
    for (b=a+1; b<A->nseq; b++)
      {
	tot+=sarseq_pair2ev (A->seq_al[a], A->seq_al[b], 1);
	n++;
      }
  return tot;
}
int sarseq_pair2ev  ( char *s1, char *s2,int mode)
{
  int l, t1, t2, t11,a, n1, n2, s;
  if ( (l=strlen (s1))!=strlen (s2))
    {
      return -1;
    }
  if (mode==2)
    {
      t1=l/2;
      t2=l/2;
      t11=l/2;
    }
  else
    {
      for (t1=t2=t11=0,a=0; a<l; a++)
	{
	  if ( mode==1)//CORRELATED
	    {
	      n1=(s1[a]=='I')?1:0;
	      n2=(s2[a]=='I')?1:0;
	    }
	  else if ( mode==0)//ANTICORRELATED
	    {
	      n1=(s1[a]=='I')?1:0;
	      n2=(s2[a]=='O')?1:0; 
	    }
	  t1+=n1;
	  t2+=n2;
	  
	  t11+=(n1&&n2)?1:0;
	}
    }
  
  
  s=(int)(100*evaluate_sar_score1 (l, t11, t1, t2));
  return s;
}
      
double evaluate_sar_score1 ( int N, int n11, int n1msa, int n1sar)
{
  double p, p1, p2;
  int n10, n01;
  
  n10=n1msa-n11;
  n01=n1sar-n11;
  
  if ( n11==0)return 0;
  //  if ( (n10)>n11 || n01>n11)return 0;
  
  p1= M_chooses_Nlog (n1msa, N) + M_chooses_Nlog (n1sar-n11, N-n1msa) + M_chooses_Nlog (n11, n1msa);
  p2=(M_chooses_Nlog (n1msa, N)+  M_chooses_Nlog (n1sar, N));
  p=(p1-p2);

  return -p;
  
}
double evaluate_sar_score2 ( int N, int n11, int n1msa, int n1sar)
{
  
  
  return n11-((n1msa-n11)+(n1sar-n11));
  
  if ( n11<n1msa) return 0;
  else if ( n11<n1sar) return 0;
  else if ( n11==N)return 0;
  return n11;
}


int benchmark_sar( int value)
{
  static int v[1000];
  static int a;

  if (a==0)
    {
      for (a=0; a< 1000; a++)v[a]=0;
      v[2]=1; 
      v[3]=2;
      v[6]=2; 
      v[7]=1;
      v[8]=2; 
      v[9]=1;
      v[10]=1; 
      v[11]=1;
      v[12]=2; 
      v[30]=2;
      v[31]=1; 
      v[32]=2;
      v[33]=1; 
      v[34]=2;
      v[35]=1;
      v[36]=1; 
      v[37]=2;
      v[43]=2; 
      v[44]=1;
      v[45]=2; 
      v[73]=2;
      v[74]=1; 
      v[75]=1;
      v[76]=2; 
      v[80]=2;
      v[81]=1; 
      v[82]=2;
      v[83]=1; 
      v[85]=2;
      v[86]=1;
      v[87]=1;
      v[88]=2; 
      v[89]=2;
      v[90]=1; 
      v[91]=2;
      v[92]=1; 
      v[93]=2;
      v[103]=2; 
      v[104]=1;
      v[105]=1; 
      v[106]=1;
      v[107]=2; 
      v[130]=2;
      v[131]=1;
      v[132]=2;
      v[133]=1;
      v[134]=1; 
      v[135]=1;
      v[136]=2; 
      v[137]=1;
      v[138]=2;
      v[271]=2;
      v[272]=1; 
      v[273]=2;
      v[281]=2; 
      v[282]=1;
      v[283]=2; 
      v[284]=1;
      v[285]=1; 
      v[286]=1;
      v[287]=2;
      v[319]=2;
      v[320]=1; 
      v[321]=1;
      v[322]=1; 
      v[323]=1;
      v[324]=2; 
      v[325]=1;
      v[326]=2; 
      v[327]=1;
      v[328]=2; 
      v[356]=2;
      v[357]=1; 
      v[358]=1;
      v[359]=2; 
      v[377]=2;
      v[378]=1;
      v[379]=2; 
      v[386]=3;
      v[388]=2;
      v[389]=1; 
      v[390]=1;
      v[391]=1; 
      v[392]=2;
      v[393]=2; 
      v[394]=2;
      v[395]=1; 
      v[396]=1;
      v[397]=2; 
      v[399]=2;
      v[400]=1; 
      v[401]=2;
      v[414]=2;
      v[415]=1;
      v[416]=2; 
      v[420]=2;
      v[421]=1; 
      v[422]=1;
      v[423]=1; 
      v[424]=2;
      v[425]=1; 
      v[426]=2;
    }
  return v[value];
}

Alignment *weight2sar (Alignment *A, Alignment *SAR, char *weight_file, int limit)
{
  int a, b, c;
  int ***weight;
  char ***list;
  float score;
  
  weight=vcalloc (SAR->nseq, sizeof (int**));
  
  
  list=file2list (weight_file, " ");

  a=b=0;
  for (a=0; a< SAR->nseq; a++)
    {
      b=c=0;
      while (list[b])
	{
	  if ( strm (list[b][1], SAR->name[a]) && atoi (list[b][3])>0)c++;
	  b++;
	}

      weight[a]=declare_int (c+1, 3);
      fprintf ( stderr, "\n%s %d", SAR->name[a], c);
      b=c=0;
      while (list[b])
	{
	  if ( strm (list[b][1], SAR->name[a]) && atoi (list[b][3])>0)
	    {
	      weight[a][c][0]=atoi(list[b][2])-1;
	      weight[a][c][1]=list[b][5][0];
	      weight[a][c][2]=atoi (list[b][3]);
	      c++;
	    }
	  b++;
	}
      weight[a][c][0]=-1;
    }
 
  for (a=0; a<A->nseq; a++)
    {
      fprintf ( stdout, ">%s\n", A->name[a]);
      for ( b=0; b< SAR->nseq; b++)
	{
	  score=seq2weighted_sar_score(A->seq_al[a], weight[b]);
	  fprintf ( stdout, "%c", (score>limit)?'I':'O');
	}
      fprintf (stdout, "\n");
    }
  myexit (EXIT_SUCCESS);
  return A;
}
  
Alignment *display_sar ( Alignment *A, Alignment *SAR, char *compound)
{
  int a,c;
  char name[100];
  
  c=name_is_in_list ( compound, SAR->name, SAR->nseq, 100);
  if ( c==-1)return A;

  for ( a=0; a< A->nseq; a++)
    {
      sprintf (name, "%s", A->name[a]);
      sprintf ( A->name[a], "%c_%s_%s", SAR->seq_al[c][a], name,compound);
    }
  return A;
}
Alignment *aln2weighted_sar_score ( Alignment *A,Alignment *SAR, char *weight_file, char *compound)
{
  
  int a, b, c=0;
  int **weight;
  
  int score;
  char reactivity;
  char ***list;

  
  if ( SAR)
    {
      c=name_is_in_list (compound, SAR->name, SAR->nseq, 100);
    }
  
  list=file2list (weight_file, " ");
  a=b=0;
  while (list[a])
    {
      if (strm (list[a][1], compound))b++;
      a++;
    }
  weight=declare_int ( b+1, 3);
  
  
  a=b=0;
  while (list[a])
    {
      if ( !strm (list[a][1], compound) || strm ("TOTPOS", list[a][1]));
      else
	{
	  weight[b][0]=atoi(list[a][2])-1;
	  weight[b][1]=list[a][5][0];
	  weight[b][2]=atoi(list[a][3]);
	  b++;
	}
      a++;
    }
  weight[b][0]=-1;
  for ( a=0; a< A->nseq; a++)
    {
      score=seq2weighted_sar_score (A->seq_al[a], weight);
      reactivity=(!SAR || c==-1)?'U':SAR->seq_al[c][a];
      
      sprintf (A->seq_comment[a], "Compound %-15s Reactivity %c SAR_SCORE %5d", compound,reactivity, (int) score);
      
    }
  return A;
}

float seq2weighted_sar_score ( char *seq, int **weight)
{
  int a, p, r, w;
  float score=0;
  
  a=0;
  while (weight[a][0]!=-1)
    {
      p=weight[a][0];
      r=weight[a][1];
      w=weight[a][2];
      
      if ( is_gap(seq[p]));
      else if ( tolower(seq[p])==r)score+=w;
      a++;
    }
  return score;
  }

Alignment * sar2simpred (Alignment *A, Alignment *SAR, char *posfile, char *compound, int L1,int L2 )
{
  int a, b, c, c1, c2;
  int **sim, **sim_ref, npred=0;
  float n11, n10, n01, n00;
  float sn, sp; 
  
  int tot_sim=0;
  int N11=1, N01=2, N10=3, NXX=4, SIM=5;
  float ***tot;
  int i1, i2;
  
  
  n11=n10=n01=n00=0;
  tot=declare_arrayN(3,sizeof (float), 10, 6, 2);
  
  sim_ref=aln2sim_mat (A, "idmat");
  if (strm (posfile, "all"))
    sim=sim_ref;
  else
    {
      Alignment *B;
      B=copy_aln ( A,NULL);
      B=extract_aln3(B,posfile);
      
      /*if (B->len_aln==0)L1=100;
      else
	L1=((B->len_aln-1)*100)/B->len_aln;
      
      if (L1<=0)L1=100;
      */
      sim=aln2sim_mat (B, "idmat");
    }
  
  for (a=0; a< A->nseq-1; a++)
    {
      for ( b=a+1; b< A->nseq; b++)
	{
	  for ( c=0; c<SAR->nseq; c++)
	    {
	      if ( (strm (compound, SAR->name[c]) || strm ( compound, "all")))
		{
		  /*if ( sim_ref[a][b]<30 || sim_ref[a][b]>60)continue;*/
		  i1=0; /*sim_ref[a][b]/10;if (i1==10)i1--;*/
		  
		  i2=sim[a][b];

		  
		  c1=(SAR->seq_al[c][a]=='I')?1:0;
		  c2=(SAR->seq_al[c][b]=='I')?1:0;
		  
		  n11=(c1 && c2)?1:0;
		  n01=(!c1 && c2)?1:0;
		  n10=(c1 && !c2)?1:0;
		  n00=(!c1 && !c2)?1:0;
		  
		  tot[i1][N11][0]+=n11;
		  tot[i1][N01][0]+=n01;
		  tot[i1][N10][0]+=n10;
		  /*tot[i1][N00][0]+=n00;*/
		  tot[i1][NXX][0]++;
		  tot[i1][SIM][0]+=sim_ref[a][b];
		  
		  if ( i2>=L1)
		    {
		      tot[i1][N11][1]+=n11;
		      tot[i1][N01][1]+=n01;
		      tot[i1][N10][1]+=n10;
		      /*tot[i1][N00][1]+=n00;*/
		      tot[i1][NXX][1]++;
		      tot[i1][SIM][1]+=sim_ref[a][b];
		    }
		}
	    }
	}
    }
  
  for (a=0; a<1; a++)
    {
      sp=(tot[a][N11][0])/(tot[a][N11][0]+tot[a][N10][0]);
      fprintf ( stdout, "\n%15s N11 %5d SP %.2f ",compound, (int)tot[a][N11][0],sp);
      sp=((tot[a][N11][1]+tot[a][N10][1])==0)?1:(tot[a][N11][1])/(tot[a][N11][1]+tot[a][N10][1]);
      sn=(tot[a][N11][0]==0)?1:(tot[a][N11][1]/tot[a][N11][0]);
      fprintf ( stdout, " N11 %5d SP %.2f SN %.2f SIM %.2f", (int)tot[a][N11][1], sp,sn, (tot[a][SIM][1]/tot[a][NXX][1]));
    }
  
  myexit (0);
  sp=((n11+n01)==0)?1:n11/(n11+n01);
  sn=((n11+n01)==0)?1:n11/(n11+n10);
  
  fprintf ( stdout, "\nLimit: %d NPRED %d AVGSIM %d SN %.2f   SP %.2f TP %d FP %d FN %d",L1, npred, tot_sim, sn, sp, (int)n11, (int)n01, (int)n10);
  myexit (EXIT_SUCCESS);
  return A;
}

Alignment * sar2simpred2 (Alignment *A, Alignment *SAR, char *seqlist, char *posfile, char *compound, int L )
{
  int a,b, c,c1, c2, p, s;
  float n11, n10, n01, n00, n, sn2, prediction,sp, n1, n0, t, entropy, Delta;
  int *rlist, *tlist, *pred, *npred, tsim, psim;
  int **sim, **sim_ref;
  int nr=0;
  int nrs;
  char *out;
  int delta_max;
  Alignment *B;
  int printall=1;

  out=vcalloc (A->nseq+1, sizeof (char));
  rlist=vcalloc ( A->nseq, sizeof (int));
  tlist=vcalloc ( A->nseq, sizeof (int));
  pred=vcalloc(2, sizeof (int));
  npred=vcalloc(2, sizeof (int));
  
  nrs=0;
  if ( strm (seqlist, "first"))
    {
      for ( a=0; a<SAR->nseq; a++)
	{
	  if ( strm ( compound, SAR->name[a]))
	    {
	      for ( b=0; b<A->nseq; b++)
		{
		  if ( SAR->seq_al[a][b]=='I')
		    {
		      fprintf ( stderr, "COMP: %s REF SEQ: %s\n", A->name[b], compound);
		      rlist[nrs]=b;
		      tlist[rlist[nrs]]=1;
		      nrs++;
		      break;
		    }
		}
	    }
	}
    }
  else if (strm (seqlist, "all"))
    {
      for ( a=0; a< A->nseq; a++)
	{
	  rlist[nrs]=a;
	  tlist[rlist[a]]=1;
	  nrs++;
	}
    }
  else if ((a=name_is_in_list ( seqlist, A->name, A->nseq, 100))!=-1)
    {
      rlist[nrs]=a;
      tlist[rlist[nrs]]=1;
      nrs++;
    }
  else
    {
      Alignment *R;
      R=main_read_aln (seqlist, NULL);
      for (a=0; a<R->nseq; a++)
	{
	  rlist[a]=name_is_in_list( R->name[a], A->name, A->nseq, 100);
	  tlist[rlist[a]]=1;
	}
      free_aln (R);
    }
  
  c=name_is_in_list ( compound, SAR->name, SAR->nseq, 100);
  
  sim_ref=aln2sim_mat (A, "idmat");
  if (strm (posfile, "all"))
    {
      sim=sim_ref;
      B=A;
    }
  else
    {
      B=copy_aln ( A,NULL);
      B=extract_aln3(B,posfile);
      sim=aln2sim_mat (B, "idmat");
    }
  
  n11=n10=n01=n00=n=n1=n0=0;
  delta_max=0;
  for (a=0; a<A->nseq; a++)
    {
      if ( tlist[a] && !strm (seqlist, "all"))
	out[a]=(SAR->seq_al[c][a]=='I')?'Z':'z';/*SAR->seq_al[c][a];*/
      else
	{
	  
	  pred[0]=pred[1]=0;
	  npred[0]=npred[1]=1;
	  c1=(SAR->seq_al[c][a]=='I')?1:0; 
	  for (nr=0,tsim=0,psim=0,b=0; b<nrs; b++)
	    {
	      if ( SAR->seq_al[c][rlist[b]]=='o');
	      else
		{
		  c2=(SAR->seq_al[c][rlist[b]]=='I')?1:0;
		  nr+=c2;
		  s=sim[a][rlist[b]];
		  tsim+=sim_ref[a][rlist[b]];
		  psim+=sim[a][rlist[b]];
		  if (s>=L)
		    {
		      pred[c2]+=s;
		      npred[c2]++;
		    }
		}
	    }
	  
	  if (c1==0)n0++;
	  else n1++;
	  t++;
	  
	  
	  Delta=pred[1]-pred[0];
	  
	  if (Delta<-delta_max){p=0;out[a]= (c1==0)?'O':'o';}
	  else if (Delta>delta_max){p=1;out[a]=(c1==1)?'I':'i';}
	  else {p=-1; out[a]=(c1==1)?'U':'u';}
	  
	  if ( p==-1);
	  else if (  p &&  c1)n11++;
	  else if (  p && !c1)n10++;
	  else if ( !p && !c1)n00++;
	  else if ( !p &&  c1)n01++;

	  if (p!=-1)n++;
	  if (printall)fprintf ( stdout, ">%-15s %d %c OVERALL_SIM:%d POSITION_SIM %d\n%s\n", B->name[a], c1, out[a],tsim/nrs,psim/nrs,B->seq_al[a]); 
	}
    }
  sp=((n11+n10)==0)?1:n11/(n11+n10);
  sn2=((n1)==0)?1:n11/n1;
  prediction=(n11+n00)/(n1+n0);
  entropy=(float)(M_chooses_Nlog (nr, nrs)/M_chooses_Nlog(nrs/2, nrs));
  
  fprintf ( stdout, ">%-15s Sp %.2f  Sn %.2f Pred %.2f E %.2f\n", compound,sp, sn2,prediction,entropy ); 
  fprintf ( stdout, "%s\n", out);
  
  myexit (EXIT_SUCCESS);
  return A;
}
/************************************************************************************/
/*                ALIGNMENT ANALYZE     : SAR FOR OR                                */
/************************************************************************************/

void display_or_help();
void display_or_help()
{
  fprintf ( stdout, "\nor_sar options:");
  fprintf ( stdout, "\n_ORCL_: Command_line in a file");
  
  fprintf ( stdout, "\n_ROTATE_  : rotate the sar matrix (if each entry is a compound rather than a sequence)");
  fprintf ( stdout, "\n_JNIT_    : number cycles of Jacknife");
  fprintf ( stdout, "\n_JNSEQ_   : Number of sequences picked up in alignment [0 to keep them all, 10 by default]");
  fprintf ( stdout, "\n_JSAR_    : Number of compounds picked up in the SAR matrix [0 to keep them all => default");
  fprintf ( stdout, "\n_JRSAR_   : Randomization of the SAR file between each JNIT iteration: S, C, R, SC, SR... (S: seq, C, column, R: residue");
  fprintf ( stdout, "\n_JRALN_   : Randomization of the ALN file between each JNIT iteration: S, C, R, SC, SR... (S: seq, C, column, R: residue");
  
  
  fprintf ( stdout, "\n_NPOS_    : Number of positions used to make the prediction [4 by default]");
  fprintf ( stdout, "\n_DEPTH_   : Depth of the motif degenerated alphabet   [Default=2]");
  fprintf ( stdout, "\n_POSFILE_ : Predefined list of positions in a file<P: pos score");
  fprintf ( stdout, "\n_RSAR_    : Randomization of the SAR file: S, C, R, SC, SR... (S: seq, C, column, R: residue");
  fprintf ( stdout, "\n_RALN_    : Randomization of the SAR file: S, C, R, SC, SR... (S: seq, C, column, R: residue");

  fprintf ( stdout, "\n_SARTHR_  : E-value Threshold for the filtered displays in comploo mode (def: 3)");

  

  fprintf ( stdout, "\n_MODE_    : jack, loo (leave one out), pos, sim, predict, self");
  fprintf ( stdout, "\n_COMPLIST_: Compound_list provided in a FASTA file. Each compound MUST correspond to the corresponding column in SAR");
  fprintf ( stdout, "\n_OUTALN_  : print the alignment corresponding to the list of positions");
  fprintf ( stdout, "\n_OUTTREE  : Make the tree (default nj) corrwesponding to _OUTALN_");
  fprintf ( stdout, "\n\n");
  exit (EXIT_SUCCESS);
}
Alignment * simple_sar_or(Alignment *inA, Alignment *inS, char *param);
void sar2jack (Alignment *A, Alignment *S, int nseq, int sarlen);


void sar2jack (Alignment *A, Alignment *S, int nseq, int sarlen)
{
  A=aln2jacknife (A, nseq, 0);
  S=reorder_aln  (S, A->name, A->nseq);
  S=aln2jacknife (S, 0, sarlen);
}
Alignment *or_scan (Alignment *A,Alignment *S, char *pmode)
{
  int l, a,ax,cx, b;
  char mode[100];
  int start, offset,w;
  int nl, *poslist;
  
 
  fprintf ( stdout, "\nPARAMETERS: %s\n", pmode);
  fprintf ( stderr, "\nPARAMETERS: %s\n", pmode);

  strget_param (pmode, "_SW_", "15", "%d",&w);
  strget_param (pmode, "_SMINW_", "5", "%d",&start);
  strget_param (pmode, "_SSTART_", "5", "%d",&offset);
  strget_param (pmode, "_SMODE_", "single", "%s",&mode);
  
  
  l=intlen (A->len_aln);
  poslist=vcalloc ( A->len_aln, sizeof (int));
  nl=0;
  
  for ( a=0; a< A->len_aln; a++)poslist[nl++]=a+1;
  
  if ( strm (mode, "single"))
    {
      for ( a=0; a<A->len_aln-2; a++)
	{
	  int c, gap, score;
	  for (gap=0,c=0; c<A->nseq; c++)gap+=(A->seq_al[c][a]=='-');
	  if ( !gap)
	    {
	      Alignment *B;
	      B=extract_aln (A, a+1, a+2);
	      B=or_sar (B, S, pmode, NO_PRINT);
	      score=B->score_aln;
	      free_aln (B);
	    }
	  else
	    {
	      score=0;
	    }
	  fprintf ( stdout, "P: %*d  S: %4d\n",l,a+1, score);
	  //fprintf ( stderr, "P: %*d  S: %4d\n",l,a+1, score);
		      
	}
    }
  else if  strm (mode, "scan")
    {
	  
      for ( ax=0; ax<nl; ax++)
	{
	  int best_score=0;
	  int best_pos=0, best_w=0, best_start=0, best_end=0;
	  a=poslist[ax];
	  for (b=start; b<=w; b++)
	    {
	      Alignment *B;
	      int pstart, pend;
	      pstart=a-b;
	      pend=a+b;
	      
	      if (pstart<1 || pstart>A->len_aln)continue;
	      if (pend<1 || pend>A->len_aln)continue;
	      B=extract_aln (A, a-w, a+w);
	      
	      B=or_sar (B, S,pmode, NO_PRINT);
	      
	      
	      if (B->score_aln>=best_score){best_score=B->score_aln; best_pos=a;best_w=b;best_start=pstart; best_end=pend;}
	      free_aln (B);
	    }
	  fprintf ( stdout, "P: %*d  I: %*d %*d Score: %5d L: %2d\n", l,best_pos, l,best_start+offset, l,best_end, best_score,(best_w*2)+1 );
	}
    }
  else if ( strm ( mode, "scan_comp"))
    {
      Alignment *NS=NULL;
      int *tresults;
      int s, n=0, p1;
      int  nbest;
      int *poscache;
      float sc, eval;
      int **index;
      int **poslist;
      int len=1;
      float best_sc, best_eval;
      int best_pos;
      char *best_word;
      int sth;
      int sev;
      
      index=aln2pos_simple (A, A->nseq);
      strget_param (pmode, "_NPOS_", "1", "%d", &len);
      strget_param (pmode, "_STH_", "0", "%d", &sth);
      strget_param (pmode, "_SEV_", "0", "%d", &sev);

      if (sev==0)for (a=0, sev=1; a<len; a++)sev*=A->len_aln;
      
      tresults=vcalloc ( A->len_aln, sizeof (int));
      poscache=vcalloc ( A->len_aln, sizeof (int));
      poslist=generate_array_int_list (len, 0, A->len_aln-1,1, &n, NULL);
      for (s=0; s<S->len_aln; s++)
	{

	  char *word;
	  NS=aln2block (S, s+1, s+2, NS);
	  fprintf ( stderr, "\nProcess: %s ...\n", get_compound_name(s, mode));
	  
	  for (best_sc=0,best_pos=0,best_word=NULL,a=0; a<n; a++)
	    {
	      for (b=0; b<len; b++)poscache[poslist[a][b]]=1;
	      word=or_id_evaluate2 (A,NS, pmode,poscache, NO_PRINT, &sc);
	      //if (word)HERE ("%s", word);
	      eval=0;
	      for (b=0; b<len; b++)poscache[poslist[a][b]]=0;
	      if (word && sc>best_sc)
		{
		  best_sc=sc;
		  best_eval=eval;
		  best_word=word;
		  best_pos=a;
		  nbest=1;
		}
	      else if ( word && sc>=best_sc)
		{
		  nbest++;
		}
	      else
		vfree (word);
	    }
	  nbest/=len;
	  for (a=0; a<=A->nseq; a++)
	    {
	

	      
	      fprintf (stdout, "\n>%s PPPP: S: %s SC: %d EV: %d NBest: %d W: %s",  get_compound_name(s, mode),(a==A->nseq)?"cons":A->name[a],(int)best_sc, (int) best_eval, nbest,best_word);
	      for ( b=0; b<len; b++)
		{
		  if ( a==A->nseq) 
		    {
		      p1=poslist[best_pos][b];
		      if (best_sc>sth && nbest<sev)tresults[p1]++;
		      p1++;
		    }
		  else 
		    {
		      p1=index[a][poslist[best_pos][b]];
		      if (p1>0)p1++;
		    }
		  fprintf ( stdout, " %d ", p1);
		  
		}
	    }
	  
	}
      for ( p1=0; p1<A->len_aln;p1++)
	fprintf ( stdout, "\n>TOT_P: %d %d", p1+1, tresults[p1]);
    }
  
  else if ( strm ( mode, "scan_comp_old"))
    {
      Alignment *NS=NULL, *BLOCK=NULL;
      int **results, *tresults;
      int s, n, p1, p2;
      int npos=1;
      int *poscache;
      
      results=declare_int (A->len_aln*A->len_aln, 3);
      tresults=vcalloc ( A->len_aln, sizeof (int));
      poscache=vcalloc ( A->len_aln, sizeof (int));
      for (s=0; s<S->len_aln; s++)
	{
	  int count;
	  NS=aln2block (S, s+1, s+2, NS);
	  fprintf ( stderr, "\nProcess: %s ...", get_compound_name(s, mode));
	  for (n=0,p1=0; p1<A->len_aln-w; p1++, count ++)
	    {
	      if ( count == 50){fprintf ( stderr, "*");count=0;}
	      for ( p2=p1; p2<A->len_aln-w; p2++, n++)
		{
		  poscache[p1]=1;
		  poscache[p2]=1;
		  
		  BLOCK=alnpos2block (A,poscache,BLOCK);
		 		
		  if ( aln2ngap(BLOCK)<=0)BLOCK=or_sar (BLOCK,NS, pmode, NO_PRINT);
		  else BLOCK->score_aln=0;
		  
		  //if ( BLOCK->score_aln>0)HERE ("P: %d %d %d", p1, p2,BLOCK->score_aln);
		  results[n][0]=p1;
		  results[n][1]=p2;
		  results[n][2]=BLOCK->score_aln;
		  poscache[p1]=poscache[p2]=0;
	    
		}
	    }
	  sort_int_inv (results, 3, 2, 0, n-1);
	  for (p1=0; p1<npos; p1++)
	    {
	      fprintf ( stdout, "\n>%s PPPP: %d %d SC: %d",  get_compound_name(s, mode), results[p1][0]+1,results[p1][1]+1, results[p1][2]);
	      fprintf ( stderr, "\n>%s PPPP: %d %d SC: %d",  get_compound_name(s, mode), results[p1][0]+1,results[p1][1]+1, results[p1][2]);
	      tresults[results[p1][0]]++;
	      tresults[results[p1][1]]++;
	    }
	}
      for ( p1=0; p1<A->len_aln;p1++)
	if ( tresults[p1])fprintf ( stdout, "\n>TOT_P: %d %d", p1+1, tresults[p1]);
    }
  
  exit (EXIT_SUCCESS);
}


ORP * or_sar_compound(Alignment *A, Alignment *S, char *mode, int print);
ORP* set_orp_name ( ORP* P, char *name);
ORP* set_orp_offset ( ORP* P,int offset);
Alignment * or_sar(Alignment *inA, Alignment *inS, char *mode, int print)
{
  char rsar[4],raln[4], sarmode[100];
  int rotate, a, b, start, end;
  ORP **R, *ORS;
  Alignment *A, *S;
  
  if ( mode && (strstr (mode, "help") || strstr (mode, "HELP")))
   display_or_help();
  
  strget_param (mode, "_SARMODE_", "single", "%s", &sarmode);//single | all
  strget_param (mode, "_RSAR_", "NO", "%s", rsar);
  strget_param (mode, "_RALN_", "NO", "%s", raln);
  strget_param (mode, "_ROTATE_", "0", "%d", &rotate);
  
  strget_param (mode, "_START_", "1", "%d", &start);start--;
  strget_param (mode, "_END_", "0", "%d", &end);
  
  A=copy_aln (inA, NULL);
  S=copy_aln (inS, NULL);
  
  if ( end==0)end=A->len_aln;
  if ( start!=0 || end!=A->len_aln)
    {
      int c;
      if ( start>=A->len_aln || end<=start || end >A->len_aln)
	printf_exit (EXIT_FAILURE, stderr, "ERROR: _START_%d and _END_%d are incompatible with the aln [len=%d][FATAL:%s]", start, end, A->len_aln, PROGRAM);
      for (b=0; b<A->nseq; b++)
	for (c=0,a=start; a<end; a++, c++)
	  A->seq_al[b][c]=A->seq_al[b][a];
      A->len_aln=c;
      for (b=0; b<A->nseq; b++)A->seq_al[b][c]='\0';
    }
  
  
  S=aln2random_aln (S, rsar);
  A=aln2random_aln (A, raln);

 

  if (rotate)
    {
      Alignment *rS;

      if (S->len_aln!=A->nseq)
	printf_exit ( EXIT_FAILURE,stderr, "ERROR: Alignment and SAR matrix are incompatible [FATAL:%s]", PROGRAM);
      
      rS=rotate_aln (S, NULL);
      
      for (a=0; a< A->nseq; a++)sprintf (rS->name[a], "%s", A->name[a]);
      free_aln (S);
      S=rS;
    }
  
  R=vcalloc ( S->len_aln+2, sizeof (ORP*));
  if (strm (sarmode, "all"))R[0]=or_sar_compound (A, S, mode, print);
  else if ( strm (sarmode, "single"))
    {
      for ( a=0; a<S->len_aln; a++)
	{
	  Alignment *NS=NULL;
	  
	  NS=aln2block (S, a+1, a+2, NS);
	  
	  R[a]=or_sar_compound (A, NS, mode, NO_PRINT);
	  set_orp_name ((R[a]), get_compound_name (a, mode));
	  set_orp_offset(R[a], start);
	  display_or_summary (R[a], mode, stdout, print);
	  
	}
    }

  
  
  ORS=combine_n_predictions (R, A, S);
  
  display_or_summary (ORS, mode, stdout, print);
  inA->score_aln=(int)(ORS->best*(float)1000);
  
  free_orp_list (R);
  free_orp(ORS);
  
  return inA;
 }
ORP* set_orp_offset ( ORP* P,int offset)
{
  if (!P) return NULL;
  else
    {
      P->offset=offset;
      return set_orp_offset(P->PR, offset);
    }
}
  
ORP* set_orp_name ( ORP* P, char *name)
{
  if (!P) return NULL;
  else
    {
      sprintf ( P->name, "%s", name);
      return set_orp_name(P->PR, name);
    }
}

ORP * combine_n_predictions (ORP**R, Alignment *A, Alignment *S)
{
  int a=0;
  ORP*N=NULL;
  while (R[a])
    {
      N=combine_2_predictions (R[a++], N, A, S);
    }
  sprintf ( N->name, "ALL");
  sprintf ( N->mode, "COMBINED");
 
  return N;
}

ORP *combine_2_predictions ( ORP*IN, ORP *TO,Alignment *A, Alignment *S)
{
  int a;
  
  if ( !TO)
    {
      TO=declare_or_prediction (IN->ncomp, IN->nseq, IN->len);
      TO->A=A;
      TO->S=S;
      TO->P=copy_aln(S, NULL);
      TO->offset=IN->offset;
      TO->ncomp=0;
    }
  
  for (a=0; a< IN->len; a++)
    {
      TO->pos[a]+=IN->pos[a];
    }
  
  TO->fp+=IN->fp;
  TO->fn+=IN->fn;
  TO->tp+=IN->tp;
  TO->tn+=IN->tn;
  rates2sensitivity (TO->tp, TO->tn, TO->fp, TO->fn, &(TO->sp), &(TO->sn), &(TO->sen2), &(TO->best));

 
  for (a=0; a<(TO->A)->nseq; a++)
    {
      (TO->P)->seq_al[a][TO->ncomp]=(IN->P)->seq_al[a][0];
      //(TO->S)->seq_al[a][TO->ncomp]=(IN->S)->seq_al[a][0];
    }
  TO->ncomp++;
  (TO->P)->len_aln=TO->ncomp;
  (TO->A)->score_aln=TO->best;
 
  return TO;
}
      
  
  
ORP * display_or_summary (ORP *CP, char *mode, FILE *fp, int print)
{
  int a;
  char *pred;
  char *exp;
  char *motif;
  Alignment *A, *P, *S;
 
  
  
  A=CP->A;
  P=CP->P;
  S=CP->S;


  
  pred=vcalloc ( P->nseq*P->len_aln*2, sizeof (char));
  exp=vcalloc ( P->nseq*P->len_aln*2, sizeof (char));
  motif=vcalloc (CP->len+1, sizeof (char));
  

  
  if (P && S)
    {
      for ( a=0; a<P->nseq; a++)
	{
	  strcat (pred,P->seq_al[a]);
	  strcat (exp, S->seq_al[a]);
	}
      CP->evalue=profile2evalue(pred, exp);
    }
  a=0;
  while ( CP->motif && CP->motif[a] && CP->motif[0][a][0])strcat (motif, CP->motif[0][a++]);
  
  if ( print==PRINT)
    {
      fprintf (fp, "\n>%-10s Mode: %s Accuracy: %6.2f E-value: %6.2f Motif: ",CP->name,CP->mode, CP->best, CP->evalue);
      
      if (motif[0])
	{
	  fprintf (fp, " %s",motif);
	}
      for ( a=0; a<CP->len; a++)
	{
	  if ( CP->pos[a]) fprintf ( fp, "\n>%-10s Mode: %s P: cons %3d SC: %4d", CP->name, CP->mode, a+1+CP->offset, CP->pos[a]);
	}
      fprintf ( fp, "\n");
    }
  vfree (pred); vfree(motif); vfree(exp);
  if (CP->PR)display_or_summary (CP->PR, mode, fp, print);
  return CP;
}
       

ORP * or_sar_compound(Alignment *A, Alignment *S, char *mode, int print)
{
  char rmode[100];
  Alignment *P=NULL;
  ORP *PR=NULL;
  
  strget_param (mode, "_MODE_", "predict", "%s", rmode);
  
  
  
  if (strm (rmode, "jack"))P=or_jack (A, S, mode);
  else if (strm (rmode, "loo")) PR=or_loo (A, S, mode,NULL, print);
  else if (strm (rmode, "comploo")) P=or_comp_loo (A, S, mode,NULL, print);
  else if ( strm (rmode, "comppos")){or_comp_pos ( A, S, mode,print);exit(0); return NULL;}
  else if ( strm (rmode, "pos"))P=or_aln2pos_aln ( A, S, mode);
  else if ( strm (rmode, "predict"))P=or_predict ( A, S, mode);
  else if ( strm (rmode, "self"))PR=or_self_predict ( A, S, mode, NULL, print);
  else if ( strm (rmode, "sim"))P=or_sim ( A, S, mode);
  
  else if ( strm (rmode, "test"))P=or_test ( A, S, mode);
  else
    {
       printf_exit (EXIT_FAILURE, stderr, "ERROR: %s is an unknown mode of or_sar  [FATAL:%s]\n", rmode,PROGRAM);
       return NULL;
    }

  if (!PR)
    {
      PR=vcalloc (1, sizeof (ORP));
      PR->P=P;
    }
  return PR;
}
Alignment * or_test ( Alignment *inA, Alignment *inS, char *mode)
{
  return inA;
}



float or_id_evaluate  ( Alignment *A, Alignment *S, char *mode, int *pos, int print)
{
  char *w;
  float score;
 
  w=or_id_evaluate2(A,S, mode, pos,print, &score);
 
  vfree (w);
  return score;
}
char* or_id_evaluate2  ( Alignment *A, Alignment *S, char *mode, int *pos, int print, float *rscore)
{
  static char **words;
  static int *plist;
  char *bword;
  
  int res, p,nl, w, c, s, exp, pred;
  int tp, tn, fp, fn;
  float sn, sp, sen2, best, best_score;
  
  if (!A)
    {free_char (words, -1);
      vfree (plist); 
      return NULL;
    }
  rscore[0]=0;
  
  plist=pos2list (pos, A->len_aln, &nl);
  words=declare_char (A->nseq, nl+1);
  bword=vcalloc (nl+1, sizeof (char));
  for (p=0; p<nl; p++)
    {
      for (s=0; s< A->nseq; s++)
	{
	  res=A->seq_al[s][plist[p]];
	  if (res=='-'){or_id_evaluate2 (NULL, NULL, NULL, NULL, 0, NULL);vfree (bword);return 0;}
	  words[s][p]=res;
	}
    }
  
  for (best_score=0,w=0; w<A->nseq; w++)
    {
      tp=fp=fn=tn=0;
      
      for (c=0; c<S->len_aln; c++)
	{
	  for (s=0; s<A->nseq; s++)
	    {
	      exp=(S->seq_al[s][c]=='I')?1:0;
	      pred=strm (words[w], words[s]);
	      if      ( exp &&  pred)tp++;
	      else if ( exp && !pred)fn++;
	      else if (!exp && !pred)tn++;
	      else if (!exp &&  pred)fp++;
	    }
	}
      rates2sensitivity (tp, tn, fp,fn,&sp, &sn, &sen2, &best);
      if ( best>best_score)
	{
	  best_score=best;
	  sprintf (bword, "%s", words[w]);
	}
    }
  rscore[0]=(float)1000*best_score;
  or_id_evaluate2 (NULL, NULL, NULL, NULL, 0, NULL);
  return bword;
}

float or_loo_evaluate2 ( Alignment *A, Alignment *S, char *mode, int *pos, int print)
{
  int c, s, p, res, sar;

  int *plist, nl;
  int tp, tn, fp, fn;
  float sn, sp, sen2, best;
  char **words, **positive, **negative;
  
  tp=tn=fp=fn=0;
  plist=pos2list (pos, A->len_aln, &nl);

  words=declare_char (A->nseq, nl+1);
  
  for (p=0; p<nl; p++)
    {
      for (s=0; s< A->nseq; s++)
	{
	  res=A->seq_al[s][plist[p]];
	  if (res=='-'){vfree (plist);free_char (words, -1); return 0;}
	  words[s][p]=res;
	}
    }
  positive=vcalloc ( A->nseq, sizeof (char*));
  negative=vcalloc ( A->nseq, sizeof (char*));
  for (c=0; c<S->len_aln; c++)
    {
      //Fill the match matrix
      for (p=0; p<nl; p++)
	{
	  for (s=0; s< A->nseq; s++)
	    {
	      sar=S->seq_al[s][c];
	      if (sar=='I')positive[s]=words[s];
	      else if ( sar=='O')negative[s]=words[s];
	    }
	}
      
      //Evaluate the scores
      for (s=0; s< A->nseq; s++)
	{
	  int pos=0, neg=0, pred;
	  sar=S->seq_al[s][c];
	  positive[s]=negative[s]=NULL;

	  if ( name_is_in_list (words[s], positive, A->nseq, nl+1)!=-1)
	    pos=1;
	  if ( name_is_in_list (words[s], negative, A->nseq, nl+1)!=-1)
	    neg=1;
	  
	  if (pos & !neg) pred=1;
	  else pred=0;
	  
	  if      ( pred  && sar=='I')tp++;
	  else if (!pred  && sar=='I')fn++;
	  else if (!pred  && sar=='O')tn++;
	  else if ( pred  && sar=='O')fp++;
	
	  if ( sar=='I')positive[s]=words [s];
	  else negative[s]=words[s];
	}
    }
  
  vfree (negative); vfree (positive);
  vfree (plist); free_char (words, -1);
  rates2sensitivity (tp, tn, fp,fn,&sp, &sn, &sen2, &best);
  
  return (float)1000*best;
}
float or_loo_evaluate ( Alignment *A, Alignment *S, char *mode, int *pos, int print)
{
  int c, s, p, res, sar;
  int **matP,**matN;
  int *plist, nl;
  int tp, tn, fp, fn;
  float sn, sp, sen2, best;
  
  tp=tn=fp=fn=0;
  plist=pos2list (pos, A->len_aln, &nl);
  matP=declare_int (nl, 256);
  matN=declare_int (nl, 256);
  
  for (c=0; c<S->len_aln; c++)
    {
      //Fill the match matrix
      for (p=0; p<nl; p++)
	{
	  for (s=0; s< A->nseq; s++)
	    {
	      res=A->seq_al[s][plist[p]];
	      sar=S->seq_al[s][c];
	      if (res=='-'){vfree (plist); free_int (matP, -1);free_int (matN, -1); return 0;}
	      if (sar=='I')matP[p][res]++;
	      if (sar=='O')matN[p][res]++;
	    }
	}
      
      //Evaluate the scores
      for (s=0; s< A->nseq; s++)
	{
	  int scoreP, scoreN;
	  int pred, valP, valN;
	  
	  sar=S->seq_al[s][c];
	  for (scoreN=0,scoreP=0,p=0; p<nl; p++)
	    {
	      res=A->seq_al[s][plist[p]];
	      
	      valP=matP[p][res]-(sar=='I')?1:0;
	      scoreP+=(valP>0)?1:0;

	      valN=matN[p][res]-(sar=='O')?1:0;
	      scoreN+=(valN>0)?1:0;
	    }

	  if ( scoreP==nl && scoreN<nl)pred=1;
	  else pred=0;
	  
	  
	  if      ( pred  && sar=='I')tp++;
	  else if (!pred  && sar=='I')fn++;
	  else if (!pred  && sar=='O')tn++;
	  else if ( pred  && sar=='O')fp++;
	}

      //reset the matrix
      for (p=0; p<nl; p++)
	{
	  for (s=0; s< A->nseq; s++)
	    {
	      res=A->seq_al[s][plist[p]];
	      sar=S->seq_al[s][c];
	      if (sar=='I')matP[p][res]=0;
	      else matN[p][res]=0;
	    }
	}
    }
  
  vfree (plist); free_int (matP, -1);free_int (matN, -1);
  rates2sensitivity (tp, tn, fp,fn,&sp, &sn, &sen2, &best);
  
  return (float)1000*best;
}
int* or_comp_pos ( Alignment *inA, Alignment *inS, char *mode,int print)
{
  Alignment *A=NULL, *S=NULL, *inS2=NULL;
  int a, b, c;
  int *main_pos, *pos=NULL;
  
  set_sar (inA, inS, mode);
  main_pos=vcalloc ( inA->len_aln, sizeof (int));
  
  inS2=copy_aln (inS, NULL);
  inS2->len_aln=1;
  
 
  
  //Run every SAR, one at a time
  for ( c=0; c< inS->len_aln; c++)
    {
      int max, p;

      fprintf ( stdout, ">%d\n", c);
      for (a=0; a< inS->nseq; a++)
	{
	  inS2->seq_al[a][0]=inS->seq_al[a][c];
	  inS2->seq_al[a][1]='\0';
	}
      
      vfree (pos);
      free_aln (S);
      free_aln (A);
      
      pos=vcalloc (inA->len_aln, sizeof (int));
      A=copy_aln (inA, NULL);
      S=copy_aln (inS2, NULL);
      set_sar (A,S, mode);
      pos=aln2predictive_positions (A, S, mode,PRINT);
      
      for (max=0,b=0; b<A->len_aln; b++)
	{
	  main_pos[b]+=pos[b];
	  if (main_pos[b]>max)
	    {
	    max=main_pos[b];
	    p=b+1;
	  }
	}

      
      for (a=0; a<A->nseq; a++)
	{
	  fprintf ( stdout, "\t");
	  for ( b=0; b<A->len_aln; b++)
	    if ( pos[b]) fprintf ( stdout, "%c", A->seq_al[a][b]);
	  fprintf ( stdout, " %c\n", inS2->seq_al[a][0]);
	}
      fprintf ( stdout, "\n\tBest: %d %d\n", p, max);
    
    }
  
  if (print==PRINT)
    {
      for ( a=0; a<inA->len_aln; a++)fprintf ( stdout, "\nP2: cons %4d %4d [FINAL]", a+1, main_pos[a]);
    }
  return main_pos;
}				
Alignment * or_comp_loo ( Alignment *inA, Alignment *inS, char *mode, int *pos,int print)
{
  int a, b,c, n;
  char **keep, **remove;
  Alignment *A, *S, *P, *P1, *SEQ, *inS2;
  int **main_pos, *compound_pos, **comp_list;
  int pos_exists=0;
  char *comp_pred, *comp_exp;
  int sar_threshold;
  
  strget_param (mode, "_SARTHRES_", "3", "%d", &sar_threshold);

  if (pos)pos_exists=1;
  set_sar (inA, inS, mode);
  P=copy_aln (inS, NULL);
  keep=declare_char (inA->nseq, MAXNAMES);
  remove=declare_char (inA->nseq, MAXNAMES);
  
  main_pos=declare_int ( inA->len_aln,4);
  comp_list=declare_int (inA->len_aln, sizeof (int*));
  inS2=copy_aln (inS, NULL);
  inS2->len_aln=1;
  
  
  comp_pred=vcalloc ( inA->nseq+1, sizeof (int));
  comp_exp=vcalloc ( inA->nseq+1, sizeof (int));
  compound_pos=NULL;
  //Run every SAR, one at a time
  for ( c=0; c< inS->len_aln; c++)
    {
      for (a=0; a< inS->nseq; a++)
	{
	  inS2->seq_al[a][0]=inS->seq_al[a][c];
	  inS2->seq_al[a][1]='\0';
	}
      vfree (compound_pos);
      compound_pos=vcalloc (inA->len_aln, sizeof (int));
      for (a=0; a<inA->nseq; a++)
	{
	  char ***motifs;
	  
	  A=copy_aln (inA, NULL);
	  S=copy_aln (inS2, NULL);
	  
	  for (n=0,b=0; b<A->nseq; b++)
	    {
	      if (b!=a)sprintf (keep[n++ ], "%s", A->name [b]);
	    }
	  sprintf ( remove[0], "%s", A->name[a]);
	  reorder_aln (A,keep, A->nseq-1);
	  
	  set_sar (A,S, mode);
	  if (!pos_exists)
	    {
	      pos=aln2predictive_positions (A, S, mode,NO_PRINT);
	    }
	  for (b=0; b<A->len_aln; b++)
	    {
	      compound_pos[b]+=pos[b];
	    }
	  
	  motifs=compounds2motifs (A, S, pos,0, mode, NO_PRINT);
	  
	  SEQ=copy_aln (inA, NULL);
	  SEQ=reorder_aln (SEQ, remove, 1);
	  
	  P1=aln2prediction (SEQ, motifs, pos);
	  comp_pred[a]=P1->seq_al[0][0];
	  comp_exp[a]=inS2->seq_al[a][0];
	  P->seq_al[a][c]=P1->seq_al[0][0];
	  
	  
	  free_aln (SEQ);
	  free_aln (S);
	  free_aln (A);
	  free_aln (P1);
	  free_arrayN( (void *)motifs, 3);
	  if (!pos_exists)vfree (pos);
	}
      if (print==PRINT)
	fprintf ( stdout, ">%-15s SC: %.2f E; %.2f\n%s\n%s\n", get_compound_name(c, mode),profile2sensitivity (comp_pred, comp_exp, NULL, NULL, NULL, NULL),profile2evalue(comp_pred, comp_exp),comp_pred, comp_exp);
      for (b=0; b<A->len_aln; b++)
	    {
	      main_pos[b][2]+=compound_pos[b];
	      if (compound_pos[b])main_pos[b][3]++;
	      if (profile2evalue(comp_pred, comp_exp)>sar_threshold)
		{
		  main_pos[b][0]+=compound_pos[b];
		  if (compound_pos[b])
		    {
		      main_pos[b][1]++;
		      comp_list[b][0]++;
		      comp_list[b]=vrealloc (comp_list[b], sizeof (int)*(comp_list[b][0]+1));
		      comp_list[b][comp_list[b][0]]=c;
		    }
		}
	    }

    }
 
  P->score_aln=(int)((float)1000*evaluate_prediction (P, inS, mode,print));
  
  if (print==PRINT)
    {
      for ( a=0; a<inA->len_aln; a++)fprintf ( stdout, "\nP: cons %4d RS: %4d RC: %5d FC: %4d %4d[FINAL]", a+1, main_pos[a][2], main_pos[a][3], main_pos[a][0], main_pos[a][1]);
      
      for ( a=0; a<inA->len_aln; a++)
	{
	  fprintf ( stdout, "\nP: cons %4d RS: %4d RC: %5d FC: %4d %4d CLIST: ", a+1, main_pos[a][2], main_pos[a][3], main_pos[a][0], main_pos[a][1]);
	  for ( c=1; c<=comp_list[a][0]; c++)
	    {
	      fprintf ( stdout, "%s ",  get_compound_name(comp_list[a][c], mode));
	    }
	  fprintf ( stdout, " [COMP_LIST]");
	}
    }
  free_int (main_pos, -1);
  free_int (comp_list, -1);
  return P;
}

ORP* or_loo ( Alignment *inA, Alignment *inS, char *mode, int *pos,int print)
{
  int a, b,  n;
  char **keep, **remove;
  Alignment *A, *S, *P, *P1, *SEQ;
  
  int pos_exists=0;
  ORP *PR;



  if (pos)pos_exists=1;
  set_sar (inA, inS, mode);
  PR=declare_or_prediction (inS->nseq, inA->nseq, inA->len_aln);
  sprintf (PR->mode, "loo ");
  P=copy_aln (inS, NULL);


  keep=declare_char (inA->nseq, MAXNAMES);
  remove=declare_char (inA->nseq, MAXNAMES);
  
    
  PR->A=inA;
  PR->P=P;
  PR->S=inS;
  
  for (a=0; a<inA->nseq; a++)
    {
      char ***motifs;
    
      A=copy_aln (inA, NULL);
      S=copy_aln (inS, NULL);
      
      for (n=0,b=0; b<A->nseq; b++)
	{
	  if (b!=a)sprintf (keep[n++ ], "%s", A->name [b]);
	}
      sprintf ( remove[0], "%s", A->name[a]);
      reorder_aln (A,keep, A->nseq-1);
      
      set_sar (A,S, mode);
      
      if (!pos_exists)
	{
	  pos=aln2predictive_positions (A, S, mode,print);
	  
	}
      
      for (b=0; b<A->len_aln; b++)
	{
	  PR->pos[b]+=pos[b];
	}
      
      motifs=compounds2motifs (A, S, pos,0, mode, print);
      
      SEQ=copy_aln (inA, NULL);
      SEQ=reorder_aln (SEQ, remove, 1);
      SEQ->nseq=1;
      
      P1=aln2prediction (SEQ, motifs, pos);
      
      
      if (print==PRINT)
	{
	  fprintf ( stdout, "\n%s\nPred: %s\nReal: %s\n", P1->name[0], P1->seq_al[0], inS->seq_al[a]);
	}
      sprintf ( P->seq_al[a], "%s", P1->seq_al[0]);
      free_aln (P1);
      
      
      free_aln (SEQ);
      free_aln (S);
      free_aln (A);

      free_arrayN( (void *)motifs, 3);
      if (!pos_exists)vfree (pos);
      
    }
  free_char (keep, -1);
  free_char (remove, -1);
  

  PR=new_evaluate_prediction (PR, mode,print);
  
    
  PR->PR=or_self_predict(inA, inS, mode,NULL, print); 
  
  if (print==PRINT)for ( a=0; a<inA->len_aln; a++)fprintf ( stdout, "\nP: cons %d %d [FINAL]", a+1, PR->pos[a]);
  
  

  return PR;
}


      
Alignment * or_jack(Alignment *inA, Alignment *inS, char *mode)
{
  int a,b;
  int n_cycles=100;
  int subnseq=10;
  int subsar=0;
  Alignment *A, *S;
  int *main_pos,*pos;
  char jrsar[10], jraln[10];
  
  strget_param (mode, "_JNIT_", "100", "%d", &n_cycles);
  strget_param (mode, "_JNSEQ_", "10", "%d", &subnseq);
  strget_param (mode, "_JNAR_", "0", "%d", &subsar);
  
  strget_param (mode, "_JRSAR_", "NO", "%s", jrsar);
  strget_param (mode, "_JRALN_", "NO", "%s", jraln);
  

 
  main_pos=vcalloc ( inA->len_aln, sizeof (int));
  for (a=0; a<n_cycles; a++)
    {
      
      A=copy_aln (inA, NULL);
      S=copy_aln (inS, NULL);
      
      S=aln2random_aln (S, jrsar);
      A=aln2random_aln (A, jraln);
      
      set_sar (A,S, mode); 
      sar2jack (A, S,subnseq,subsar);
      
      

      pos=aln2predictive_positions (A, S,mode, PRINT);
      
      for (b=0; b<A->len_aln; b++)main_pos[b]+=pos[b];
      vfree (pos);

    }
  display_pos (A, S, main_pos, mode);
  
 
  return inA;
}

Alignment * display_pos (Alignment *A, Alignment *S, int *pos,char *mode)
{
  Alignment *B;
  int a, b;
  int **index;

  int intl;

  intl=intlen (A->len_aln);
  index=aln2pos_simple (A, A->nseq);
  B=copy_aln (A,NULL);
  B->len_aln=0;
  for ( a=0; a<A->len_aln; a++)
    fprintf ( stdout, "\nP: cons %*d %*d S: %4d [DISPLAY_FULL_POS]", intl,a+1,intl, a+2, pos[a]);
  fprintf ( stdout, "\n\n");
  for (a=0; a<A->len_aln; a++)
    {
      if (pos[a])
	{
	  for ( b=0; b<A->nseq; b++)
	    {
	      B->seq_al[b][B->len_aln]=A->seq_al[b][a];
	      if (index[b][a]>0)fprintf ( stdout, "\nP: %s %d %d S: %d [DISPLAY_POS]",A->name[b], index[b][a], index[b][a]+1, pos[a]);
	    }
	  B->len_aln++;
	  fprintf ( stdout, "\nP: cons %d %d S: %d [DISPLAY_POS]", a+1, a+2, pos[a]);
	}
    }
  fprintf ( stdout, "\n");
  for (a=0; a<B->nseq; a++)B->seq_al[a][B->len_aln]='\0';
  return B;
}
Alignment * or_aln2pos_aln (Alignment *A, Alignment *S, char *mode)
{
  Alignment *B;

  int *pos;
  char outaln[100], outtree[100];
  

  strget_param (mode, "_OUTALN_", "NO", "%s", outaln);
  strget_param (mode, "_OUTTREE_", "NO", "%s", outtree);
 
  set_sar (A, S, mode);
  pos=aln2predictive_positions (A, S,mode, PRINT);
  
  B=display_pos (A, S, pos, mode);
  
  
  if (!strm(outaln, "NO")) vfclose (output_aln (B, vfopen (outaln, "w")));
  if (!strm(outtree, "NO"))vfclose (print_tree (aln2tree(B), "newick", vfopen (outtree, "w")));
 
  return B;
}
Alignment * or_sim(Alignment *A, Alignment *S, char *mode)
{
  //Predict all the sequences that are not both in inS and inA
  int *pos;
  
  set_sar (A, S, mode);
  pos=aln2predictive_positions (A, S,mode, PRINT);
  fprintf ( stdout, "R: %.3f", pos2sim (A,S, pos));
  
  exit (EXIT_SUCCESS);
  return A;
}
ORP* or_self_predict(Alignment *A, Alignment *S, char *mode,int *pos, int print)
{
  //Predict all the sequences that are not both in inS and inA
  Alignment *P;
  char ***motifs;
  

  int a;

  int pre_set_pos=0;
  ORP *PR;

  
  set_sar (A, S, mode);
  PR=declare_or_prediction (S->nseq, A->nseq, A->len_aln);
  sprintf (PR->mode, "self");
  PR->A=A;
  PR->S=S;
  
  if (!pos)
    {
      pos=aln2predictive_positions (A, S,mode,print);
      pre_set_pos=0;
    }
  else
    pre_set_pos=1;

  for (a=0; a< A->len_aln; a++)
    PR->pos[a]=pos[a];
  

  PR->motif=motifs=compounds2motifs (A, S, pos,0, mode, print);
  P=PR->P=aln2prediction (A, motifs, pos);
  
  if (!pre_set_pos)vfree (pos);
  
  PR=new_evaluate_prediction (PR, mode,print);
  return PR;
}


Alignment * or_predict(Alignment *inA, Alignment *inS, char *mode)
{
  //Predict all the sequences that are not both in inS and inA
  Alignment *P, *A, *S, *T;
  char ***motifs;
  int *pos;

  int a, b;


  
  
  

  A=copy_aln (inA, NULL);
  S=copy_aln (inS, NULL);
  set_sar (A, S, mode);
  
  pos=aln2predictive_positions (A, S,mode,PRINT);
  motifs=compounds2motifs (A, S, pos,0, mode, PRINT);
  T=get_prediction_target (inA, inS, mode);
  

  P=aln2prediction (T, motifs, pos);
  //recall=evaluate_prediction (S, P, mode);
  for ( a=0; a<P->len_aln; a++)
    {
      for (b=0; b<P->nseq; b++)
	{
	  if (tolower(P->seq_al[b][a])=='i')fprintf (stdout, "\n>%20s %20s %c", T->name [0],get_compound_name (a, mode), P->seq_al[b][a]);
	}
    }
  fprintf ( stdout, "\n");
  return P;
}

Alignment *get_prediction_target (Alignment *A, Alignment *S, char *param)
{
  char **name;
  int n, a;
  Alignment *T;
  
  T=copy_aln (A, NULL);
  name=declare_char (A->nseq, 100);
  for (n=0,a=0; a< A->nseq; a++)
    {
      if ( name_is_in_list (A->name[a], S->name, S->nseq, 100)==-1)
	{
	  sprintf (name[n++], "%s", A->name[a]);
	}
    }
  T=reorder_aln (T,name, n);
  return T;
}

Alignment *set_sar (Alignment *A, Alignment *S, char *param)
{
  char **name;
  int n, a;
  
  name=declare_char (A->nseq, 100);
  for (n=0,a=0; a< A->nseq; a++)
    {
      if ( name_is_in_list (A->name[a], S->name, S->nseq, 100)!=-1)
	{
	  sprintf (name[n++], "%s", A->name[a]);
	}
    }
  A=reorder_aln (A,name, n);
  S=reorder_aln (S,name, n);
  free_char (name, -1);
  return S;
}  

ORP* new_evaluate_prediction  (ORP *PR, char *mode, int print)
{
  int a,b, i, r, p;
  int tp, tn, fp, fn;
  float sn, sp, sen2, best;
  float tot_best_seq=0;
  float tot_best_comp=0;
  Alignment *P, *R;
  
  int ns=0;
  float *recall;

  
  P=PR->P;
  R=PR->S;
  
  recall=vcalloc (P->len_aln, sizeof (float));
  if (P->len_aln!=R->len_aln)
    {
      HERE ("Mismatch between number of compounds in prediction and reference");
      exit (0);
    }
  if (print==PRINT)fprintf ( stdout, "\n");
  
  for (a=0; a<P->nseq; a++)
    {
      tp=tn=fp=fn=0;
      if ((i=name_is_in_list (P->name[a], R->name, R->nseq, 100))!=-1)
	{
	  
	  for (b=0;b<P->len_aln; b++)
	    {
	      r=R->seq_al[i][b];
	      p=P->seq_al[a][b];
	      
	      if ( p=='I' && r=='I')tp++;
	      else if ( p=='I' && r=='O')fp++;
	      else if ( p=='O' && r=='I')fn++;
	      else if ( p=='O' && r=='O')tn++;
	    }
	  rates2sensitivity (tp, tn, fp, fn, &sp, &sn, &sen2, &best);
	  if (print==PRINT)fprintf (stdout, ">%-s sp: %.2f sn: %.2f sn2: %.2f best: %.2f [SEQ]\n",P->name[a], sp, sn, sen2, best);
	  if ( best>0)
	    {
	      ns++;
	      tot_best_seq+=best;
	    }
	}
    }
  if (ns)
    {
      tot_best_seq/=ns;
    }
  if (print==PRINT)fprintf ( stdout, ">TotSeq sp: %.2f N: %d[SEQ]\n",tot_best_seq, ns);
  
  tot_best_comp=0;
  for (ns=0,b=0; b<P->len_aln; b++)
    {
      tp=tn=fp=fn=0;
      for (a=0; a<P->nseq;a++) 
	{
	  if ((i=name_is_in_list (P->name[a], R->name, R->nseq, 100))!=-1)
	    {
	      r=R->seq_al[i][b];
	      p=P->seq_al[a][b];
	      
	      if ( p=='I' && r=='I'){PR->tp++;tp++;}
	      else if ( p=='I' && r=='O'){PR->fp++;fp++;}
	      else if ( p=='O' && r=='I'){PR->fn++;fn++;}
	      else if ( p=='O' && r=='O'){PR->tn++;tn++;}
	    }
	}
      rates2sensitivity (tp, tn, fp, fn, &sp, &sn, &sen2, &best);
      
      if (print==PRINT) fprintf (stdout, ">%-25s sp: %.2f sn: %.2f sen2: %.2f best: %.2f [COMP]\n",get_compound_name (b, mode), PR->sp, PR->sn, PR->sen2,PR->best);
      if ( best>0)
	{
	  ns++;
	  tot_best_comp+=best;
	}
    }
  
  if (ns)
    {
      tot_best_comp/=ns;
    }
  rates2sensitivity (PR->tp, PR->tn, PR->fp,PR->fn,&(PR->sp), &(PR->sn), &(PR->sen2), &(PR->best));
  if (print==PRINT)fprintf ( stdout, ">FullTot sp: %.2f sn: %.2f sen2: %.2f best: %.2f  N: %d[COMP]\n", PR->sp, PR->sn, PR->sen2,PR->best, ns);
  P->score_aln=(int)((float)1000*(PR->best));
  return PR;
}
float  evaluate_prediction  (Alignment *R, Alignment *P, char *mode, int print)
{
  int a,b, i, r, p;
  int tp, tn, fp, fn;
  int tot_tp, tot_tn, tot_fp, tot_fn;
  float sn, sp, sen2, best;
  float tot_sp=0;
  float tot_sn=0;
  float tot_sen2=0;
  float tot_best_seq=0;
  float tot_best_comp=0;
  float tot_best=0;
  
  int ns=0;
  float *recall;
 
  
 
  
  recall=vcalloc (P->len_aln, sizeof (float));
  if (P->len_aln!=R->len_aln)
    {
      HERE ("Mismatch between number of compounds in prediction and reference");
      exit (0);
    }
  if (print==PRINT)fprintf ( stdout, "\n");
  for (a=0; a<P->nseq; a++)
    {
      tp=tn=fp=fn=0;
      if ((i=name_is_in_list (P->name[a], R->name, R->nseq, 100))!=-1)
	{
	  
	  for (b=0;b<P->len_aln; b++)
	    {
	      r=R->seq_al[i][b];
	      p=P->seq_al[a][b];
	      
	      if ( p=='I' && r=='I')tp++;
	      else if ( p=='I' && r=='O')fp++;
	      else if ( p=='O' && r=='I')fn++;
	      else if ( p=='O' && r=='O')tn++;
	    }
	  rates2sensitivity (tp, tn, fp, fn, &sp, &sn, &sen2, &best);
	  if (print==PRINT)fprintf (stdout, ">%-s sp: %.2f sn: %.2f sn2: %.2f best: %.2f [SEQ]\n",P->name[a], sp, sn, sen2, best);
	  if ( best>0)
	    {
	      ns++;
	      tot_best_seq+=best;
	      tot_sn+=sn;
	      tot_sp+=sp;
	      tot_sen2+=sen2;
	    }
	}
    }
  if (ns)
    {
      tot_best_seq/=ns;
      tot_sn/=ns;
      tot_sp/=ns;
      tot_sen2/=ns;
    }
  if (print==PRINT)fprintf ( stdout, ">Tot sp: %.2f sn: %.2f sen2: %.2f best: %.2f  N: %d[SEQ]\n", tot_sp, tot_sn, tot_sen2,tot_best_seq, ns);
  
  tot_fp=tot_fn=tot_tp=tot_tn=0;
  tot_sp=tot_sn=tot_sen2=tot_best_comp=0;
  for (ns=0,b=0; b<P->len_aln; b++)
    {
      tp=tn=fp=fn=0;
      for (a=0; a<P->nseq;a++) 
	{
	  if ((i=name_is_in_list (P->name[a], R->name, R->nseq, 100))!=-1)
	    {
	      r=R->seq_al[i][b];
	      p=P->seq_al[a][b];
	      
	      if ( p=='I' && r=='I'){tot_tp++;tp++;}
	      else if ( p=='I' && r=='O'){tot_fp++;fp++;}
	      else if ( p=='O' && r=='I'){tot_fn++;fn++;}
	      else if ( p=='O' && r=='O'){tot_tn++;tn++;}
	    }
	}
      rates2sensitivity (tp, tn, fp, fn, &sp, &sn, &sen2, &best);
      
      if (print==PRINT) fprintf (stdout, ">%-25s sp: %.2f sn: %.2f sen2: %.2f best: %.2f [COMP]\n",get_compound_name (b, mode), sp, sn, sen2,best);
      recall[b]=sen2;
      if ( best>0)
	{
	  ns++;
	  tot_best_comp+=best;
	  tot_sn+=sn;
	  tot_sp+=sp;
	  tot_sen2+=sen2;
	}
    }
  
  if (ns)
    {
      tot_best_comp/=ns;
      tot_sn/=ns;
      tot_sp/=ns;
      tot_sen2/=ns;
    }
  rates2sensitivity (tot_tp, tot_tn, tot_fp,tot_fn,&tot_sp, &tot_sn, &tot_sen2, &tot_best);
  if (print==PRINT)fprintf ( stdout, ">FullTot sp: %.2f sn: %.2f sen2: %.2f best: %.2f  N: %d[COMP]\n", tot_sp, tot_sn, tot_sen2,tot_best, ns);
  return tot_best;
}
	  


Alignment * aln2prediction (Alignment *A,char ***motif, int *pos)
{
  int a, b,nc, nl;
  int *list;
  char **array, **sar;
  Alignment *R;
  Sequence *S;
  nc=read_array_size ((void *)motif, sizeof (char***));
  
  
  list=pos2list (pos, A->len_aln, &nl);

  
  array=declare_char (A->nseq, nl+1);
  sar=declare_char(A->nseq, nc+1);
  for (a=0; a<A->nseq; a++)
    {
      for (b=0; b<nl; b++)
	array[a][b]=A->seq_al[a][list[b]];
    }

  for (a=0; a<nc; a++)
    {
      for (b=0; b<A->nseq; b++)
	{

	  sar[b][a]=(match_motif (array[b], motif[a]))?'I':'O';
	}
    }

  
  S=fill_sequence_struc (A->nseq,sar,A->name);
  R=seq2aln (S, NULL, KEEP_GAP);
  free_sequence (S, S->nseq);
  free_char (sar, -1);
  vfree (list);
  free_char (array, -1);
  return R;
}

int *   file2pos_list (Alignment *A, char *posfile)
{
  char ***file;
  int **index;
  int *pos;
  int i, n, p;

  //pos_file:<P: seqname pos score\n>
  //          1  2       3   4
  
  
  if ( !check_file_exists (posfile))
    {
      printf_exit ( EXIT_FAILURE, stderr, "ERROR: Could not read posfile %s\n", posfile);
    }
  
  file=file2list (posfile, " ");
 
  index=aln2inv_pos (A);
  pos=vcalloc ( A->len_aln, sizeof (int));
  
  n=0;
  while (file[n])
    {
      
      if ( !strm (file[n][1], "P:"));
      else
	{
	  if ( (strm (file[n][2], "cons")))
	    p=atoi(file[n][3])-1;
	  else 
	    {
	      i=name_is_in_list ( file[n][2], A->name, A->nseq, MAXNAMES+1);
	      if (i!=-1)
		p=index[i][atoi(file[n][3])]-1;
	      else p=-1;
	    }
	  if (p!=-1)pos[p]+=atoi(file[n][4]);
	}
      n++;
    }

  
  free_int (index, -1);
  free_arrayN ( (char **)file, 3);
  return pos;
} 
int *   aln2predictive_positions (Alignment *A, Alignment *B, char *mode, int print)
{
  char posmode[100];
  
  if (!mode) return NULL;
  HERE ("%s", mode);
  strget_param (mode, "_POSMODE_", "scan", "%s", posmode);
  if ( strm (posmode, "mat"))return aln2predictive_positions_mat (A, B, mode, print);
  else if ( strm (posmode, "scan")) return aln2predictive_positions_scan (A, B, mode, print);
  else
    {
      printf_exit (EXIT_FAILURE,stderr, "ERROR: %s is an unknown _POSMODE_ mode",posmode);
      return NULL;
    }
}

int *   aln2predictive_positions_mat (Alignment *A, Alignment *B, char *mode, int print)
  {
    int a, b, c,gap,  res1, res2, sar1, sar2, npos, s, idscore;
    float id1,id2,id3,nid1,nid2,nid3;
    int **pos, *fpos;
    pos=declare_int (A->len_aln,2);
    fpos=vcalloc ( A->len_aln, sizeof (int));
    
    strget_param (mode, "_NPOS_", "2", "%d", &npos);
    for ( a=0; a< A->len_aln; a++)
      {
	pos[a][0]=a;
	id1=id2=id3=nid1=nid2=nid3=0;
	for ( gap=0,b=0; b<A->nseq; b++)gap+=(A->seq_al[b][a]=='-');
	if ( gap>0){pos[a][1]=0;continue;}
	
	for (s=0; s<B->len_aln; s++)
	  {
	    for ( gap=0,b=0; b<A->nseq-1; b++)
	      {
		sar1=B->seq_al[b][s];
		res1=A->seq_al[b][a];
		
		for ( c=b+1; c<A->nseq; c++)
		  {
		    sar2=B->seq_al[c][s];
		    res2=A->seq_al[c][a];
		    
		    idscore=(res1==res2)?1:0;
		    if ( sar1 == 'I' && sar2=='I'){id1+=idscore;nid1++;}
		    else if ( sar1 =='0' && sar2=='0'){id2+=idscore;nid2++;}
		    else {id3+=idscore; nid3++;}
		   
		  }
	      }	
	    id1=(nid1==0)?1:id1/nid1;
	    id2=(nid1==0)?1:id2/nid2;
	    id3=(nid3==0)?1:id3/nid3;
	    pos[a][1]=(int)((float)1000*id1*(1-id3));
	    
	  }
      }
    
    sort_int (pos, 2,1, 0, A->len_aln-1);
    for ( a=MAX(0,(A->len_aln-npos));a<A->len_aln; a++)
      {
	fpos[pos[a][0]]=1;
      }
    
    free_int (pos, -1);
    return fpos;
  }
int *   aln2predictive_positions_scan (Alignment *A, Alignment *B, char *mode, int print)
{
  int a, b, c, best_pos,nl, nplist=0, max, posw;
  float best_score, score;
  static int *list, *tpos,**plist,*array;
  int *pos;


  char posfile[100];
  char predmode[100];
  char target_posfile[100];

  

  if (!A)
    {
      vfree (list);
      vfree (tpos);
      
      free_int (plist, -1);
      vfree (array);
      return NULL;
    }

  strget_param (mode, "_PREDMODE_", "ID", "%s", predmode);
  strget_param (mode, "_POSW_", "1", "%d", &posw);
  strget_param (mode, "_NPOS_", "2", "%d", &max);
  strget_param (mode, "_POSFILE_", "NO", "%s", posfile);
  strget_param (mode, "_TPOSFILE_", "NO", "%s", target_posfile);
  
  if ( !strm(posfile, "NO"))return file2pos_list (A,posfile);
  if ( !strm(target_posfile, "NO"))tpos=file2pos_list (A,target_posfile);
  else
    {
      tpos=vcalloc (A->len_aln, sizeof (int));
      for (a=0; a<A->len_aln; a++)tpos[a]=1;
    }

  //Declare the positions that are going to be scanned
  

  if (posw==1)
    {
      plist=declare_int (A->len_aln, 2);
      nplist=0;
      for (a=0; a<A->len_aln; a++)
	{
	  if(tpos[a])
	    {
	      plist[nplist][0]=1;
	      plist[nplist][1]=a;
	      nplist++;
	    }
	}
    }
  else if ( posw==2)
    {
      nplist=0;
      plist=declare_int (A->len_aln*A->len_aln, 3);
      for (a=0; a<A->len_aln; a++)
	for (b=0; b<A->len_aln; b++)
	  {
	    plist[nplist][1]=a;
	    plist[nplist][2]=b;
	    plist[nplist][0]=2;
	    nplist++;
	  }
    }
  else if ( posw==3)
    {
      nplist=0;
      plist=declare_int (A->len_aln*A->len_aln*A->len_aln, 3);
      for (a=0; a<A->len_aln; a++)
	for (b=0; b<A->len_aln; b++)
	  {
	    plist[nplist][1]=a;
	    plist[nplist][2]=b;
	    plist[nplist][3]=c;
	    
	    
	    plist[nplist][0]=3;
	    nplist++;
	  }
    }

  
  pos=vcalloc ( A->len_aln, sizeof (int));
  if (max==0)max=A->len_aln;
  else if ( max==-1)
    {
      for (a=0; a<A->len_aln; a++)if (tpos[a]){pos[a]=1;}
      aln2predictive_positions_scan (NULL, NULL, NULL, 0);
      return pos;
    }

  

  pos=vcalloc ( A->len_aln, sizeof (int));
  list=vcalloc (A->len_aln, sizeof (int));
  nl=0;
  
  

  for (a=0; a< max; a++)
    {
      int previous_best_pos=-1;
      for (best_score=-9999,best_pos=0,b=0; b<nplist; b++)
	{
	  for (c=0; c<nl; c++)pos[list[c]]=1;
	  for (c=1; c<=plist[b][0]; c++)pos[plist[b][c]]=1;
	
	  if (strm(predmode, "R"))score=sar_aln2r(A,B,pos,0);
	  else if ( strm (predmode, "ID"))
	    {
	      score  =or_id_evaluate (A, B, mode, pos,NO_PRINT);
	    }
	  else if ( strm (predmode, "BP2"))score =or_loo_evaluate2 (A, B, mode, pos,NO_PRINT);
	  else
	    {
	      HERE ("Unknown mode: %s", predmode);
	      exit (0);
	    }
	  if ( score>best_score)
	    {
	      best_score=score;
	      best_pos=b;
	    }
	  for (c=1; c<=plist[b][0]; c++)pos[plist[b][c]]=0;
	  
	}
      if (best_pos==previous_best_pos)break;
      else previous_best_pos=best_pos;
      
      //update the best_pos_list
      for (b=1; b<=plist[best_pos][0]; b++)
	list[nl++]=plist[best_pos][b];
      
      
      if ( print==PRINT)
	{
	  for (b=0; b<nl; b++) pos[list[b]]=1;
	  fprintf ( stdout, "\nP_current: ");
	  for ( c=1; c<=plist[best_pos][0]; c++)fprintf ( stdout, "%d ",plist[best_pos][c]+1);
	  fprintf ( stdout, " S: %.3f D: %d R:%d", best_score, (int)sar_aln2delta(A,B, pos,0), nl);
	}
      for (b=0; b<nl; b++) pos[list[b]]=0;    
    }

  for (a=0; a<nl; a++) 
    pos[list[a]]=1;
  if (print==PRINT)fprintf ( stdout, "\nR_best: %.3f with %d pos" ,best_score, nl);
  
  aln2predictive_positions_scan (NULL, NULL, NULL, 0);
  
  return pos;
}

char *** compounds2motifs (Alignment *A, Alignment *B, int *pos, int depth, char *mode, int print)
{
  char ***motifs;
  int a;

  motifs=vcalloc (B->len_aln, sizeof (char**));
  for (a=0; a<B->len_aln; a++)
    {

      motifs[a]=compound2motif (A, B, pos, depth, a, mode, print);
    }
  
  return motifs;
}
char ** compound2regexp_motif (Alignment *A, Alignment *B, int *pos, int depth, int c, char *mode, int print);
char ** compound2word_motif (Alignment *A, Alignment *B, int *pos, int depth, int c, char *mode, int print);

char ** compound2motif (Alignment *A, Alignment *B, int *pos, int depth, int c, char *mode, int print)
{
  char mmode[100];
 
  strget_param (mode, "_MOTIFMODE_", "word", "%s", mmode); //words, regexp
  if ( strm (mmode, "regexp"))return compound2regexp_motif (A,B,pos, depth, c, mode, print);
  else if ( strm (mmode, "word"))return compound2word_motif (A,B,pos, depth, c, mode, print);
 else return NULL;}
char ** compound2word_motif (Alignment *A, Alignment *B, int *pos, int depth, int c, char *mode, int print)
{
  int a,l;
  char *word, **motif;
  float score;
  

  word=or_id_evaluate2 (A, B, mode, pos,print, &score);
  if ( !word) return NULL;
  l=strlen (word);
  
  motif=declare_char (l+1, 2);
  for (a=0; a<l; a++)motif[a][0]=word[a];
  
  if (print==PRINT)
    {
      fprintf ( stdout, "\nMotifCompound %25s best: %.2f  motif: ", get_compound_name(c, mode),score);
      for (a=0; a<l; a++) 
	{
	  fprintf ( stdout, "[%2s]",motif[a]); 
	}
    }
  vfree (word);
  return motif;
}
    

    
char ** compound2regexp_motif (Alignment *A, Alignment *B, int *pos, int depth, int c, char *mode, int print)
{
  static Alignment *I;
  static Alignment *O;
  int a, b, o, i;

  float tp,tn,fp,fn,best, sp, sn, sen2, best_sn, best_sp, best_sen2;
  float best_pred=-1;
  int best_motif=0;



  char ***alp=NULL;
  int *alp_size=NULL;

  char *motif_file;
  int n;
  char **m, **m2;

  char *buf=NULL;
  char *best_buf;
  FILE *fpp;
  
  free_aln (I);
  free_aln (O);
  
  I=copy_aln(A, NULL);
  O=copy_aln(A, NULL);
  
  if (depth==0)
    strget_param (mode, "_DEPTH_", "2", "%d", &depth);
  
  I->nseq=O->nseq=I->len_aln=O->len_aln=0;
  for (a=0; a<A->len_aln; a++)
    {
      if (pos[a])
	{
	  for (i=o=0,b=0; b<A->nseq; b++)
	    {
	      if ( is_gap(A->seq_al[b][a]))return 0;
	      if (B->seq_al[b][c]=='I')I->seq_al[i++][I->len_aln]=A->seq_al[b][a];
	      else O->seq_al[o++][O->len_aln]=A->seq_al[b][a];
	    }
	  I->len_aln++;
	  O->len_aln++;
	}
    }
  
  if (O->len_aln==0 || I->len_aln==0) return 0;
  O->nseq=o;
  I->nseq=i;
  for (a=0; a<o; a++)O->seq_al[a][O->len_aln]='\0';
  for (a=0; a<i; a++)I->seq_al[a][I->len_aln]='\0';
  
  if (!I->nseq) return NULL;
  

 
  best_pred=best_motif=best_sn=best_sp=best_sen2=0;
 
  motif_file=vtmpnam (NULL);
  
  n=0;
  if (depth>0)
    {
       alp=vcalloc ( sizeof (char**), I->len_aln);
       alp_size= vcalloc ( I->len_aln, sizeof (int));
       for (a=0; a<I->len_aln; a++)
	 {
	   char *col;
	   alp[a]=string2alphabet ( (col=aln_column2string (I,a)),depth, &alp_size[a]);
	   vfree (col);
	 }
       generate_array_string_list (I->len_aln, alp, alp_size, &n, motif_file, OVERLAP);
    }
  else
    {
      int *used;
      char r;
      
      used=vcalloc (256, sizeof (int));
      fpp=vfopen (motif_file,"w");
      for (a=0;a<I->len_aln; a++)
	{
	  for (b=0; b<I->nseq; b++)
	    {
	      r=I->seq_al[b][a];
	      if (!used[(int)r]){fprintf (fpp, "%c", r);used[(int)r]=1;}
	    }
	  for (b=0; b<I->nseq; b++)
	    {
	      r=I->seq_al[b][a];
	      used[(int)r]=0;
	    }
	  fprintf (fpp, " ");
	}
      fprintf (fpp, "\n");
      vfree (used);
      vfclose (fpp);
      
      n=1;
      depth=I->nseq;
    }
		 
  buf=vcalloc (2*(I->len_aln*depth)+1, sizeof (char));
  best_buf=vcalloc (2*(I->len_aln*depth)+1, sizeof (char));
  fpp=vfopen (motif_file, "r");
  
  for (a=0; a<n; a++)
    {
      
      buf=vfgets (buf, fpp);
      m2=string2list (buf);
      m2++;
      
      
      tp=tn=fp=fn=0;
      
      for (b=0; b<I->nseq; b++)
	{
	  if (match_motif (I->seq_al[b], m2))tp++;
	  else fn++;
	}
      for (b=0; b<O->nseq; b++)
	{
	  if (match_motif (O->seq_al[b], m2))fp++;
	  else tn++;
	}
      rates2sensitivity (tp, tn, fp, fn, &sp, &sn, &sen2, &best);
      
      if (best>= best_pred)
	{
	  best_pred=best;
	  best_sp=sp;
	  best_sen2=sen2;
	  best_sn=sn;
	  sprintf (best_buf, "%s", buf);
	}
      m2--;
      free_char (m2, -1);
    }
  vfclose (fpp);
  if (print==PRINT)fprintf ( stdout, "\nMotifCompound %25s sp: %.2f sn: %.2f sen2: %.2f best: %.2f  motif: ", get_compound_name(c, mode), best_sp, best_sn, best_sen2, best_pred);
  m2=string2list (best_buf);
  m=declare_char (I->len_aln+1, depth+1);
  
  for (a=0; a<I->len_aln; a++) 
	{
	  sprintf (m[a], "%s", m2[a+1]);
	  if (print==PRINT) fprintf ( stdout, "[%2s]",m[a]); 
	}
  if (print==PRINT)fprintf ( stdout, " N-motifs %d", n);
  free_char (m2, -1);
  
  if (alp)free_arrayN((void ***) alp, 3);
  if (alp_size)vfree (alp_size);
  vfree (buf); vfree(best_buf);
  
  return m;
}

double pos2sim (Alignment *A, Alignment *B, int *pos)
{
  return sar_aln2r (A, B,pos, PRINT);
}
double  sar_aln2r (Alignment *A, Alignment *B, int *pos, int print)
{
  int a, b, c, d,r1, r2, n, score, sim;
  double *r, result;
  static double **slist;
  int declare=0;
  static int **M;
 

  
  if (!M)M=read_matrice ("blosum62mt");
  if (!slist)
    {
      int maxslist;
      maxslist=A->nseq*A->nseq*10;
      slist=declare_double (maxslist, 2);
    }

  if (pos==NULL)
    {

      declare=1;
      pos=vcalloc ( A->len_aln+1, sizeof (int));
      for (a=0; a<A->len_aln; a++)pos[a]=1;
      pos[a]=-1;
    
    }
  
  for (n=0,a=0; a< A->nseq-1; a++)
    {
      
      for (b=a+1; b<A->nseq; b++)
	{

	  
	  for (sim=d=0,c=0; c<A->len_aln; c++)
	    {
	      
	      if (pos[c]==0)continue;
	      
	      r1=A->seq_al[a][c];
	      r2=A->seq_al[b][c];
	      if (is_gap(r1) || is_gap(r2))return 0;
	      
	      sim+=M[r1-'A'][r2-'A']*pos[c];
	      d+=MAX((M[r1-'A'][r1-'A']),(M[r2-'A'][r2-'A']));
	    }
	  sim=(d==0)?0:(100*sim)/d;
	  score=(int)get_sar_sim(B->seq_al[a], B->seq_al[b]);
	  slist[n][0]=(double)sim;
	  slist[n][1]=(double)score;
	  if (print==PRINT)fprintf ( stdout, "SIM: %d %d [%s %s]\n", sim, score, A->name[a], A->name[b]);
	  n++;
	}
    }

  r=return_r(slist, n);
  for (a=0; a<n; a++)slist[a][0]=slist[a][1]=0;
  result=r[0];
  vfree (r);
  if (declare) vfree (pos);
  
  return result;
}


double sar_aln2delta (Alignment *A, Alignment *B, int *pos, int print)
{
  static Alignment *I;
  static Alignment *O;
  int a, b, c, o, i;
  double delta=0;
  if (!I)
    {
      I=copy_aln(A, NULL);
      O=copy_aln(A, NULL);
    }
  

  
  for (c=0; c<B->len_aln; c++)
    {
    
      I->nseq=O->nseq=I->len_aln=O->len_aln=0;
      for (a=0; a<A->len_aln; a++)
	{
	  if (pos[a])
	    {
	      for (i=o=0,b=0; b<B->nseq; b++)
		{
		  if ( is_gap(A->seq_al[b][a]))return 0;
		  if (B->seq_al[b][c]=='I')I->seq_al[i++][I->len_aln]=A->seq_al[b][a];
		  else O->seq_al[o++][O->len_aln]=A->seq_al[b][a];
		}
	      I->len_aln++;
	      O->len_aln++;
	    }
	}
      if (O->len_aln==0 || I->len_aln==0) return 0;
      O->nseq=o;
      I->nseq=i;
      for (a=0; a<o; a++)O->seq_al[a][O->len_aln]='\0';
      for (a=0; a<i; a++)I->seq_al[a][I->len_aln]='\0';

      delta+=aln2sim(I,"blosum62mt")-aln2sim(O, "blosum62mt");
      
    }

  return delta;
}

char * get_compound_name (int c, char *mode)
{
  static int isset;
  static Alignment *S;
  static char *lname;

  if (!isset)
    {
      char *comp_list;
      isset=1;
      lname=vcalloc (100, sizeof (char));
      
      if (!mode);
      else
	{
	  strget_param (mode, "_COMPLIST_", "NO", "%s", comp_list=vcalloc (100, sizeof (char)));
	  if (strm(comp_list, "NO"));
	  else
	    {
	      S=main_read_aln (comp_list, NULL);
	      vfree (comp_list);
	    }
	}
    }
  if (!S || c>=S->nseq)sprintf (lname, "%d", c);
  else
    {
      sprintf (lname, "%s", S->name [c]);
    }
  return lname;
}
ORP * declare_or_prediction ( int ncomp, int nseq, int len)
{
  ORP *P;
  P=vcalloc ( 1, sizeof (ORP));
  P->ncomp=ncomp;
  P->nseq=nseq;
  P->len=len;
  P->PR=NULL;
  
  P->pos=vcalloc (len+1, sizeof (int));
  
  return P;
}

void free_orp_list ( ORP**P)
{
  int a=0;
  while (P[a])
    {
      free_orp(P[a++]);
    }
}
void free_orp ( ORP*P)
{
  if (!P) return;
  free_aln (P->A);
  free_aln (P->S);
  free_aln (P->P);
  vfree (P->pos);
  free_arrayN((void **)P->motif, 3);
  if (P->PR)free_orp(P->PR);
  vfree (P);
}














  
/*********************************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**********************************/