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

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

static void make_p_ptrs(int *tptr, int *pl, int naseq, int l);
static void make_n_ptrs(int *tptr, int *pl, int naseq, int len);
static void put_frag(int fs, int v1, int v2, int flen);
static int frag_rel_pos(int a1, int b1, int a2, int b2);
static void des_quick_sort(int *array1, int *array2, int array_size);
static void pair_align(int seq_no, int l1, int l2);


/*
*	Prototypes
*/

/*
*	 Global variables
*/
/*extern int *seqlen_array;
  extern char **seq_array;
  extern int  dna_ktup, dna_window, dna_wind_gap, dna_signif; params for DNA 
 extern int prot_ktup,prot_window,prot_wind_gap,prot_signif; params for prots 
  extern int 	nseqs;
  extern Boolean 	dnaflag;
  extern double 	**tmat;
  extern int 	max_aa;
  extern int  max_aln_length;
*/

static int *seqlen_array;
static char **seq_array;
 
static int 	nseqs;
static int 	dnaflag;
static int      max_aln_length;
static int      max_aa;

static int 	next;
static int 	curr_frag,maxsf;
static int 	**accum;
static int 	*diag_index;
static char 	*slopes;

int ktup,window,wind_gap,signif;    		      /* Pairwise aln. params */
int *displ;
int *zza, *zzb, *zzc, *zzd;

static Boolean percent=1;


static void make_p_ptrs(int *tptr,int *pl,int naseq,int l)
{
	static int a[10];
	int i,j,limit,code,flag;
	int residue;
	
	/*tptr--> pointer to the last occurence of the same residue or ktuple:

	  abcdeabef
	  
	  tptr: 0 0 0 0 0 1 2 5 0
	  pl[a]=6
	  pl[b]=7
	*/

	
	for (i=1;i<=ktup;i++)
           a[i] = (int) pow((double)(max_aa+1),(double)(i-1));

	limit = (int)   pow((double)(max_aa+1),(double)ktup);
	for(i=1;i<=limit;++i)
		pl[i]=0;
	for(i=1;i<=l;++i)
		tptr[i]=0;
	
	

	for(i=1;i<=(l-ktup+1);++i) {
		code=0;
		flag=FALSE;
		for(j=1;j<=ktup;++j) {
			residue = seq_array[naseq][i+j-1];
			if((residue<0) || (residue > max_aa)){
				flag=TRUE;
				break;
			}
			code += ((residue) * a[j]);
		}
		if(flag)
			continue;
		++code;
		if(pl[code]!=0)tptr[i]=pl[code];
		pl[code]=i;
	}
}


static void make_n_ptrs(int *tptr,int *pl,int naseq,int len)
{
	static int pot[]={ 0, 1, 4, 16, 64, 256, 1024, 4096 };
	int i,j,limit,code,flag;
	int residue;
	
	limit = (int) pow((double)4,(double)ktup);
	
	for(i=1;i<=limit;++i)
		pl[i]=0;
	for(i=1;i<=len;++i)
		tptr[i]=0;
	
	for(i=1;i<=len-ktup+1;++i) {
		code=0;
		flag=FALSE;
		for(j=1;j<=ktup;++j) {
			residue = seq_array[naseq][i+j-1];
			if((residue<0) || (residue>4)){
				flag=TRUE;
				break;
			}
			code += ((residue) * pot[j]);  /* DES */
		}
		if(flag)
			continue;
		++code;
		if(pl[code]!=0)
			tptr[i]=pl[code];
		pl[code]=i;
	}
}


static void put_frag(int fs,int v1,int v2,int flen)
{
	int end;
	accum[0][curr_frag]=fs;
	accum[1][curr_frag]=v1;
	accum[2][curr_frag]=v2;
	accum[3][curr_frag]=flen;
	
	if(!maxsf) {
		maxsf=1;
		accum[4][curr_frag]=0;
		return;
	}
	
        if(fs >= accum[0][maxsf]) {
		accum[4][curr_frag]=maxsf;
		maxsf=curr_frag;
		return;
	}
	else {
		next=maxsf;
		while(TRUE) {
			end=next;
			next=accum[4][next];
			if(fs>=accum[0][next])
				break;
		}
		accum[4][curr_frag]=next;
		accum[4][end]=curr_frag;
	}
}


static int frag_rel_pos(int a1,int b1,int a2,int b2)
{
	int ret;
	
	ret=FALSE;
	if(a1-b1==a2-b2) {
		if(a2<a1)
			ret=TRUE;
	}
	else {
		if(a2+ktup-1<a1 && b2+ktup-1<b1)
			ret=TRUE;
	}
	return ret;
}


static void des_quick_sort(int *array1, int *array2, int array_size)
/*  */
/* Quicksort routine, adapted from chapter 4, page 115 of software tools */
/* by Kernighan and Plauger, (1986) */
/* Sort the elements of array1 and sort the */
/* elements of array2 accordingly */
/*  */
{
	int temp1, temp2;
	int p, pivlin;
	int i, j;
	int lst[50], ust[50];       /* the maximum no. of elements must be*/
								/* < log(base2) of 50 */

	lst[1] = 1;
	ust[1] = array_size;
	p = 1;

	while(p > 0) {
		if(lst[p] >= ust[p])
			p--;
		else {
			i = lst[p] - 1;
			j = ust[p];
			pivlin = array1[j];
			while(i < j) {
				for(i=i+1; array1[i] < pivlin; i++)
					;
				for(j=j-1; j > i; j--)
					if(array1[j] <= pivlin) break;
				if(i < j) {
					temp1     = array1[i];
					array1[i] = array1[j];
					array1[j] = temp1;
					
					temp2     = array2[i];
					array2[i] = array2[j];
					array2[j] = temp2;
				}
			}
			
			j = ust[p];

			temp1     = array1[i];
			array1[i] = array1[j];
			array1[j] = temp1;

			temp2     = array2[i];
			array2[i] = array2[j];
			array2[j] = temp2;

			if(i-lst[p] < ust[p] - i) {
				lst[p+1] = lst[p];
				ust[p+1] = i - 1;
				lst[p]   = i + 1;
			}
			else {
				lst[p+1] = i + 1;
				ust[p+1] = ust[p];
				ust[p]   = i - 1;
			}
			p = p + 1;
		}
	}
	return;

}





static void pair_align(int seq_no,int l1,int l2)
{
	int pot[8],i,j,l,m,flag,limit,pos,tl1,vn1,vn2,flen,osptr,fs;
	int tv1,tv2,encrypt,subt1,subt2,rmndr;
	int residue;
	
	if(dnaflag) {
		for(i=1;i<=ktup;++i)
			pot[i] = (int) pow((double)4,(double)(i-1));
		limit = (int) pow((double)4,(double)ktup);
	}
	else {
		for (i=1;i<=ktup;i++)
           		pot[i] = (int) pow((double)(max_aa+1),(double)(i-1));
		limit = (int) pow((double)(max_aa+1),(double)ktup);
	}
	
	tl1 = (l1+l2)-1;
	
	for(i=1;i<=tl1;++i) {
		slopes[i]=displ[i]=0;
		diag_index[i] = i;
	}
	

/* increment diagonal score for each k_tuple match */
/* Attempt at guessing the best band by looking at identities*/

	for(i=1;i<=limit;++i) 
	        {
		vn1=zzc[i];
		while(TRUE) 
		        {
			if(!vn1) break;
			vn2=zzd[i];
			while(vn2 != 0) 
			        {
				osptr=vn1-vn2+l2;
				++displ[osptr]; /*PLUG THE Pos Dependant Scheme Here!!!! (For Id only)*/
				vn2=zzb[vn2];
				}
			vn1=zza[vn1];
			}
		}

/* choose the top SIGNIF diagonals */

	des_quick_sort(displ, diag_index, tl1);

	j = tl1 - signif + 1;
	if(j < 1) j = 1;
 
/* flag all diagonals within WINDOW of a top diagonal */

	for(i=tl1; i>=j; i--) 
		if(displ[i] > 0) {
			pos = diag_index[i];
			l = (1  >pos-window) ? 1   : pos-window;
			m = (tl1<pos+window) ? tl1 : pos+window;
			for(; l <= m; l++) 
				slopes[l] = 1;
		}

	for(i=1; i<=tl1; i++)  displ[i] = 0; /*reset the diagonals score*/

	
	next=curr_frag=maxsf=0;	
	for(i=1;i<=(l1-ktup+1);++i) 
	        {
		encrypt=flag=0;
		for(j=1;j<=ktup;++j)
		{residue = seq_array[seq_no][i+j-1];if((residue<0) || (residue>max_aa)){flag=TRUE; break;}encrypt += ((residue)*pot[j]);}
		if(flag) continue;
		else flag=FALSE;
		
		++encrypt;	
		vn2=zzd[encrypt];
      
		/*now trying to match i-ktup and vn2-ktup*/
		while(TRUE) 
		        {
			if(!vn2) 
			        {
				flag=TRUE;
				break;
				}
			osptr=i-vn2+l2;      /*osptr=Diagonal under investigation*/
			if(slopes[osptr]!=1) /*Get the next diagonal if that one is not flagged*/
			        {
				vn2=zzb[vn2];
				continue;
			        }
			flen=0;
			fs=ktup;
			next=maxsf;	
			
			/* A-loop*/
			while(TRUE) 
			    {
			    if(!next)
				  {
				  ++curr_frag;
				  if(curr_frag>=2*max_aln_length) 
					  {

					  return;
					  }
				  displ[osptr]=curr_frag;
				  put_frag(fs,i,vn2,flen); /*sets the coordinates of the fragments*/
				  }
			    else 
				  {
				  tv1=accum[1][next];
				  tv2=accum[2][next];
				  if(frag_rel_pos(i,vn2,tv1,tv2)) 
				        {
					if(i-vn2==accum[1][next]-accum[2][next]) 
					    {
					    if(i>accum[1][next]+(ktup-1))
						fs=accum[0][next]+ktup;
					    else  
					          {
						  rmndr=i-accum[1][next];
						  fs=accum[0][next]+rmndr;
					          }
					    flen=next;
					    next=0;
					    continue;
					    }
					else 
					    {
					    if(displ[osptr]==0)
					    subt1=ktup;
					    else 
					       {
					       if(i>accum[1][displ[osptr]]+(ktup-1))
					       subt1=accum[0][displ[osptr]]+ktup;
					       else 
					          {
						  rmndr=i-accum[1][displ[osptr]];
						  subt1=accum[0][displ[osptr]]+rmndr;
						  }
					       }
					    subt2=accum[0][next]-wind_gap+ktup;
					    if(subt2>subt1) 
					        {
						flen=next;
						fs=subt2;
						}
					    else 
					        {
						flen=displ[osptr];
						fs=subt1;
						}
					    next=0;
					    continue;
					    }
				           
					}
				  else 
				        {
					next=accum[4][next];
					continue;
					}
				  }
			    break;
			}
		/*
		* End of Aloop
		*/
		
			vn2=zzb[vn2];
		}
	}

}		 

int ** show_pair(int istart, int iend, int jstart, int jend, int *in_seqlen_array, char **in_seq_array, int dna_ktup, int dna_window, int dna_wind_gap, int dna_signif,int prot_ktup, int prot_window,int prot_wind_gap,int prot_signif, int in_nseqs,int in_dnaflag, int in_max_aa, int in_max_aln_length  )
{
	int i,j,dsr;
	double calc_score;
	int **tmat;
	
	seqlen_array=vcalloc ( in_nseqs+1, sizeof(int));
	for ( i=0; i< in_nseqs; i++)seqlen_array[i+1]=in_seqlen_array[i];
	
	
	seq_array=declare_char ( in_nseqs+1, in_max_aln_length);
	for ( i=0; i< in_nseqs; i++)sprintf (seq_array[i+1], "%s",in_seq_array[i]); 

	
	nseqs=in_nseqs;
	dnaflag=in_dnaflag;
	max_aa=in_max_aa;
	max_aln_length=in_max_aln_length;

	
	tmat=declare_int ( nseqs+1, nseqs+1);
	accum=declare_int( 5, 2*max_aln_length+1);

	displ      = (int *) vcalloc( (2*max_aln_length +1), sizeof (int) );
	slopes     = (char *)vcalloc( (2*max_aln_length +1) , sizeof (char));
	diag_index = (int *) vcalloc( (2*max_aln_length +1) , sizeof (int) );

	zza = (int *)vcalloc( (max_aln_length+1),sizeof (int) );
	zzb = (int *)vcalloc( (max_aln_length+1),sizeof (int) );

	zzc = (int *)vcalloc( (max_aln_length+1), sizeof (int) );
	zzd = (int *)vcalloc( (max_aln_length+1), sizeof (int) );

        if(dnaflag) {
                ktup     = dna_ktup;
                window   = dna_window;
                signif   = dna_signif;
                wind_gap = dna_wind_gap;
        }
        else {
                ktup     = prot_ktup;
                window   = prot_window;
                signif   = prot_signif;
                wind_gap = prot_wind_gap;
        }

	for(i=istart+1;i<=iend;++i) 
	        {
		if(dnaflag)
			make_n_ptrs(zza,zzc,i,seqlen_array[i]);
		else
			make_p_ptrs(zza,zzc,i,seqlen_array[i]);
		for(j=MAX(jstart+1, i+1);j<=jend;++j) 
		        {
			    if (i!=j)
			       {
				   if(dnaflag)
				       make_n_ptrs(zzb,zzd,j,seqlen_array[j]);
				   else
				       make_p_ptrs(zzb,zzd,j,seqlen_array[j]);

				   pair_align(i,seqlen_array[i],seqlen_array[j]);

				   if(!maxsf)
				       calc_score=0.0;
				   else {
				       calc_score=(double)accum[0][maxsf];
				       if(percent) {
					   dsr=(seqlen_array[i]<seqlen_array[j]) ?
					       seqlen_array[i] : seqlen_array[j];
					   calc_score = (calc_score/(double)dsr) * 100.0;
				       }
				   }

				   tmat[i-1][j-1] = (int)(((100.0 - calc_score)/100.0)*1000);
				   tmat[j-1][i-1] = (int)(((100.0 - calc_score)/100.0)*1000);
				   fprintf ( stderr, "\r[%d %d]=> %.2f",i, j, (float)calc_score );
			       }
			}
		}

	free_int ( accum, -1);

	vfree(displ);
	vfree(slopes);
        vfree(diag_index);

	vfree(zza);
	vfree(zzb);
	vfree(zzc);
	vfree(zzd);
	return tmat;
}

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