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

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

#include "dp_lib_header.h"

int domain_pair_wise (Alignment *A,int*in_ns, int **in_l_s,Constraint_list *CL )
	{
/*******************************************************************************/
/*                SEQ_DOMAIN DP                                                    */
/*                                                                             */
/*	makes DP between the the ns[0] sequences and the ns[1]                 */
/*                                                                             */
/*	for MODE, see the function  get_dp_cost                                */
/*******************************************************************************/
	
	int scale, gop, gep, maximise; 
	int a, b, i, j,l,x;
	int best_j;

	int lenal[2], len;
	int sub;
	int match;
	
	int *ns, **l_s;


	int *f;
	int *pf;

	int *dd;
	int *dd_len;
	
	int * e;
	int *pe;
	int * e_len;
	int *pe_len;
	    
	int fop;
        int **pos0;
	
	int  **al=NULL;
	int  **pos_al=NULL;


	
	int ala,LEN;
	char *buffer;
	char *char_buf;


/*trace back variables       */
	TRACE_TYPE *buf_trace=NULL;
	TRACE_TYPE **trace;
	TRACE_TYPE k;
	TRACE_TYPE *tr;
	int **result_aln;
	int nseq;
/*Test Varaibles*/
	int score;

/*Prepare l_s and ns*/

	
	ns=vcalloc( 2, sizeof (int));
	l_s=vcalloc ( 2, sizeof (int*));
	
	ns[0]=in_ns[1];
	ns[1]=in_ns[0];
	  
	l_s[0]=in_l_s[1];
	l_s[1]=in_l_s[0];
	lenal[0]=strlen (A->seq_al[l_s[0][0]]);
	lenal[1]=strlen (A->seq_al[l_s[1][0]]);
	len=lenal[0]+lenal[1]+2;
/********************************/	
gop=(CL->gop)*SCORE_K;
gep=(CL->gep)*SCORE_K;
maximise=CL->maximise;
scale=(lenal[1]*SCORE_K*(CL->moca)->moca_scale);
/*******************************/

	

/*DO MEMORY ALLOCATION FOR DP*/
	

	


	buf_trace=vcalloc ( len, sizeof (TRACE_TYPE));	
	buffer=vcalloc ( 2*len, sizeof (char));	
	
        al    =declare_int  (2, 2*len);  
	pos_al=declare_int  (2, 2*len);
	result_aln=declare_int (1,len);
	char_buf= vcalloc (2*len, sizeof (char));	

	f     =vcalloc (len, sizeof (int));
       pf     =vcalloc (len, sizeof (int));
	e     =vcalloc (len, sizeof (int));
       pe     =vcalloc (len, sizeof (int));
	e_len =vcalloc (len, sizeof (int));
       pe_len =vcalloc (len, sizeof (int));
	dd    =vcalloc (len, sizeof (int));
	dd_len=vcalloc (len, sizeof (int));

	
	trace=declare_int (lenal[0]+2, lenal[1]+2);

/*END OF MEMORY ALLOCATION*/
	
	
		/*
		0(s)   +(dd)
		  \      |
		   \     |
		    \    |
		     \   |
		      \  |
		       \ |
		        \|
		-(e)----O
		*/ 
	
	pos0=aln2pos_simple ( A,-1, ns, l_s);	

	for ( i=0; i<=lenal[0]+1; i++)
	    {	    
	    tr=trace[i];

	    
	    for ( sub=0,j=0; j<=lenal[1]; j++)
	        {
		if      (i==0 && j==0){tr[j]=1;pe[j]=dd[j]=gop;}
		else if (i==0)        {e[j]=pe[j]=dd[j]=gop;dd_len[j]=e_len[j]=pe_len[j]=f[j]=pf[j]=0;tr[j]=-1;}
		else if (j==0)
		     {
		     for (f[j]=pf[0],best_j=0,a=1; a<=lenal[1]; a++)
		         {
			 if (f[j]!=MAX(pf[a]+scale,f[j]))
			     {
			     f[j]=pf[a]+scale;
			     best_j=a;
			     }			
			 }
		     

		     dd    [j]=e[j]=pe[j]=gop;
		     dd_len[j]=e_len[j]=0;
		     tr    [j]=best_j;
		     
		     }
		else if (i>lenal[0]);
		else
		     {
		     sub=(CL->get_dp_cost) (A, pos0, ns[0], l_s[0], i-1, pos0, ns[1], l_s[1],j-1,CL);
		     
		  
		     
		     match=pf[j-1]+sub;
		     if (a_better_than_b(pf[j]+gop, pe[j]+gep,maximise))
			 {
			 e    [j]=pf[j]+gop;
			 e_len[j]=1;
			 }
		     else 
		         {
			 e    [j]=pe    [j]+gep;
			 e_len[j]=pe_len[j]+1;
			 }
		     
		     
		     if (a_better_than_b(f[j-1]+gop, dd[j-1]+gep,maximise))
		        {
			dd    [j]=f[j-1]+gop;
			dd_len[j]=1;
			}
		     else 
		        {
			dd    [j]=dd    [j-1]+gep;
		        dd_len[j]=dd_len[j-1]+1  ;
			}
		     


		     if ( sub!=UNDEFINED)
			 {			 
			 
			 f[j] =best_int(4,maximise,&fop,e[j],match,dd[j],f[0]);
			 fop-=1;
			 
			 if      (fop==-1)fop= e_len[j]*fop;
			 else if (fop==1 )fop=dd_len[j]*fop;
			 else if (fop==2 )fop=UNDEFINED;
			 
			 
			 
			 }
		     else
		         {
			 dd[j]=e[j]=match=-10000;			 
			 f[j]=f[0];
			 fop=UNDEFINED;
			 }
		     pe    [j]  =e    [j];
		     pf    [j-1]=f    [j-1];
		     pe_len[j]  =e_len[j];
		     tr[j]      =fop;	
		     }
	
		}
	    
	    pf    [j-1]=f    [j-1];
	    pe    [j]  =e    [j];
	    pe_len[j]  =e_len[j];
	    }
	score=f[0];
        i=lenal[0]+1;
	j=0;
	ala=0;
	
	
	while (i!=0)
	      {
	     

	      k=trace[i][j];
	      if (j==0 && i<=lenal[0])
	         {
		 
		 pos_al[0][ala]=i;
		 pos_al[1][ala]=j;
		 al[0][ala]=MATCH;
		 al[1][ala]=UNALIGNED;
		 i--;
		 j=k;
		 ala++;
		 }
	      else if ( j==0 && i>lenal[0])
	         {
		 j=k;
		 i--; 
		 
		 }
	      else if (k==0)
	         {		
		 pos_al[0][ala]=i;
		 pos_al[1][ala]=j;
		 al[0][ala]=MATCH;
		 al[1][ala]=MATCH;
		 i--;
		 j--;
		
		 ala++;
		 }		
	      else if (k!=UNDEFINED && k<0 && i>0)
	         {
		 for (x=0; x< -k && i>0; x++)
		     {
		     pos_al[0][ala]=i;
		     pos_al[1][ala]=j;    
		     al[0][ala]=MATCH;
		     al[1][ala]=GAP;		 
		     i--;
		     ala++;
		     }
		 }
	      else if (k!=UNDEFINED && k>0 && j>0)
	         {
		 for ( x=0; x< k && j>0; x++)
		     {
		     pos_al[0][ala]=i;
		     pos_al[1][ala]=j;    
		     al[0][ala]=GAP;
		     al[1][ala]=MATCH;
		     j--;
		     ala++;
		     }
		 }
	      else if ( k==UNDEFINED){j=0;}
	      
	      }

	
	LEN=ala;	

	invert_list_int ( pos_al[0], LEN);
	invert_list_int ( pos_al[1], LEN);		
	invert_list_int (     al[0], LEN);
	invert_list_int (     al[1], LEN);	

	
	/*O: TARGET SEQUENCE  (long)*/
	/*1: PATTERN SEQUENCE (short)*/

	
	
	for ( b=0; b<lenal[1]; b++)result_aln[0][b]=-1;
	for (l=0, nseq=0, a=0; a< LEN; a++)
	    {
	    i=pos_al[0][a];
	    j=pos_al[1][a];
	    
	    if (al[1][a]==UNALIGNED && l>0)
	             {
		     result_aln=realloc_int ( result_aln, read_size_int ( result_aln,sizeof (int*)),len, 1, 0); 	   
		     nseq++;
		     l=0;
		     for ( b=0; b<lenal[1]; b++)result_aln[nseq][b]=-1;
		     }
	    
	    else if (al[1][a]==MATCH && al[0][a]==MATCH ){l++;result_aln[nseq][j-1]=i-1;}
	    else if (al[1][a]==MATCH && al[0][a]==GAP   ){l++;result_aln[nseq][j-1]=-1;}
	    
	    }
	if ( l>0)nseq++;




	A=domain_match_list2aln ( A,ns,l_s,result_aln,nseq,lenal[1]); 

	
	vfree (f);
	vfree (pf);
	vfree (e);
	vfree (pe);
	vfree (e_len);
	vfree (pe_len);
	vfree (dd_len);
	vfree (dd);
	free_int (pos0, -1);
	vfree (buffer);
	vfree (char_buf); 
	vfree (buf_trace);
	free_int ( pos_al, -1);
        pos_al=NULL;
	  
	free_int (     al, -1);

	free_int (trace,-1);
	free_int ( result_aln, -1);
	return score;
	}


Alignment *domain_match_list2aln ( Alignment *A,int *ns,int **l_s,int **ml, int nseq, int len)
           {

	     /*
	       function documentation: start
	       
	       This function edits the alignment given the results obtained by DP
	       ns: ns[0]->number of sequences serarched (TARGET)
	           ns[1]->number of sequences in the pattern (PATTERN SEQ)
	       l_s:   
		   l_s[0]->list of sequences in the TARGET...
	       
	       nseq: number of occurences of PATTERN in TARGET
	       len:  length of the PATTERN
	       
	       ml: detail of the nseq matches
	           ml[x][y]=k-> residue k of TARGET matches residue y of pattern
	                     -> k=-1 means a gap;
			     
	       NOTE: This implementation can only match ONE target sequence with the PATTERN
	       The Pattern can either be one sequence or a profile.

	       function documentation: end
	     */

	     
	     

	   int a, b, c, d, e;
	   Alignment *B=NULL;
	   int **new_ml;
	   int  *max_ml;
	   int  *start_ml;
	   int tot_nseq;
	   int max_len,seq;
	   char *buf;
	   
	   if ( len==0 || nseq==0)
	      {
	      A->nseq=0;
	      A->len_aln=0;
	      }
	   else
	      {
	      B=copy_aln(A, B);
	      /*1 Extract the sequence used as a pattern, put it on the top*/
	      
	    
	      A=shrink_aln (A, ns[1], l_s[1]);
	      A=realloc_aln2(A, ns[1]+ns[0]*nseq,len+A->len_aln+1);
	      
	      
	      
	      new_ml  =declare_int ( nseq, 3*len);
	      max_ml  =vcalloc ( nseq, sizeof (int));
	      for ( a=0; a<nseq; a++)max_ml[a]=-1;

	      start_ml=vcalloc ( nseq, sizeof (int));

	      for ( b=0,a=0; a< len; a++, b+=3)
	          {
		  for ( max_len=0,c=0; c< nseq; c++)
		      {
		      if ( max_ml[c]<0 && ml[c][a]<0)
		         {
			 new_ml[c][b]=ml[c][a];
		         new_ml[c][b+1]=0;
			 }
		      else if ( ml[c][a]>=0)
		         {		       
			 new_ml[c][b]=ml[c][a];
			 if ( max_ml[c]<0) start_ml[c]=max_ml[c]=ml[c][a];
			 for ( d=a+1; d<len; d++){if ( ml[c][d]>=0){ max_ml[c]= ml[c][d];break;}}
			 if (max_ml[c]!=new_ml[c][b])
			    {
			    new_ml[c][b+1]=max_ml[c]-new_ml[c][b]-1;
			    }
			 max_len=MAX( max_len, new_ml[c][b+1]);
			 }
		      else
		         {
			 new_ml[c][b]=ml[c][a];
			 new_ml[c][b+1]=0;
			 }
		      }
	       
		  for ( c=0; c< nseq; c++){new_ml[c][b+2]=max_len;}
		  }
	  
	      tot_nseq=ns[1]+ns[0]*nseq;
	   
	      for ( a=0, b=0; a< len ;a++)
	          {
		 
	          /*1: Place the Match Column*/
		  for ( c=0; c< ns[1]; c++)
		      {
		      A->seq_al[c][b]=B->seq_al[l_s[1][c]][a];
		      A->seq_al[c][b+1]='\0';		   
		      }
		  for ( e=0,c=ns[1]; c<tot_nseq; c+=ns[0],e++)
		      for ( d=0; d< ns[0]; d++)
		          {
			  if ( new_ml[e][3*a]!=-1)
			    A->seq_al[c+d][b]=B->seq_al[l_s[0][d]][new_ml[e][3*a]];
			  else 
			      A->seq_al[c+d][b]='-';
			  A->seq_al[c+d][b+1]='\0';
		       
			  }
		  b++;

	         /*2: Add the Gaps before the next_column*/
		  if ( new_ml[0][3*a+2]>0)
	             {
		     for ( c=0; c< ns[1]; c++)
		         {
			 buf=generate_null(new_ml[0][3*a+2]);
			 strcat ( A->seq_al[c],buf);
			 vfree (buf);
			 }
		     for (e=0,c=ns[1];c< tot_nseq; c+=ns[0], e++)
		         {
			   buf=extract_char (B->seq_al[l_s[0][0]], new_ml[e][3*a]+1, new_ml[e][3*a+1]);
			   strcat ( A->seq_al[c],buf);
			   vfree (buf);
			   buf=generate_null(new_ml[e][3*a+2]-new_ml[e][3*a+1]);
			   strcat ( A->seq_al[c],buf);
			   
			   vfree (buf);
			   
			 }
		     }
		  b+=new_ml[0][3*a+2];
		  }
	      
	      for (e=0,a=ns[1]; a< tot_nseq; a+=ns[0],e++)
	          {
		  for ( b=0; b<ns[0]; b++)
		      {
		      seq=l_s[0][b];
		      A->order[a+b][0]=B->order[seq][0];
		      A->order[a+b][1]=B->order[seq][1];
		      for ( c=0; c<start_ml[e];c++)A->order[a+b][1]+=!is_gap(B->seq_al[seq][c]);
		      sprintf ( A->name[a+b], "Repeat_%d", a+b);
		      }
		  }
	
	      free_aln(B);
	      A->nseq=tot_nseq;
	      A->len_aln=strlen ( A->seq_al[0]);
	      
	      }
	   return A;
	   }
Alignment * domain_seq2domain (Constraint_list *CL,int scale,int gop,int gep,Alignment *SEQ_DOMAIN, Alignment *TARGET) 
	   {	
	   static Alignment *A;
	   int *n_groups;
	   int **group_list;
	   int a,b,c;
	   
	   A=copy_aln (TARGET, A);
	   A=stack_aln( A, SEQ_DOMAIN);

	   
	   n_groups=vcalloc ( 2, sizeof (int));
	   group_list=declare_int (2, A->nseq);
	   
	   n_groups[0]=TARGET->nseq;
	   n_groups[1]=SEQ_DOMAIN->nseq;
	   for (c=0, a=0; a< 2; a++)
	       {
	       for (b=0; b< n_groups[a]; b++, c++)
	           {
		   group_list[a][b]=c;
		   }
	       }	   
	   A->score_aln=domain_pair_wise (A, n_groups, group_list,CL);
	   
	   SEQ_DOMAIN=copy_aln (A, SEQ_DOMAIN);
	   vfree (n_groups);
	   free_int (group_list,-1);
	   return SEQ_DOMAIN;
	   }
	   

/******************************COPYRIGHT NOTICE*******************************/
/*© Centro de Regulacio Genomica */
/*and */
/*Cedric Notredame */
/*Fri Feb 18 08:27:45 CET 2011 - Revision 596. */
/*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*******************************/