#ifndef ALIGN_H
#define ALIGN_H

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

#define MAX_ALPHABETSIZ 25
#define MAXSEQSIZE      120000
int ALPHABETSIZ = MAX_ALPHABETSIZ;              /* the alphabet size for compositional alignment */




/* added by Gary Benson 1/2003 */
typedef struct
{
    char*   sequence1side;
    char*   sequence2side;
	char*   matchboundaries;
    int     sequence1start;
    int     sequence1end;
    int     sequence2start;
    int     sequence2end;
    int     score;
    int     length;
}   COMPOSITIONALIGNPAIR;

typedef struct
{
	int C[MAX_ALPHABETSIZ];
} COMPOSITIONVECTOR;

/* #pragma pack is MSVC++ specific */

#pragma pack(push,1)

typedef struct
{
    int		 score;
    char     direction;
}   SD;

#pragma pack(pop)



/************************************************************************
*                       - function prototypes -
*************************************************************************/

/* added by Gary Benson 1/2003 */
COMPOSITIONALIGNPAIR* GetCompositionGlobalAlignPair(char* seq1, char* seq2, int len1,
                            int len2, int alpha, int indel, int* submatrix, int* M);

/* added by Gary Benson 1/2003 */
COMPOSITIONALIGNPAIR* GetCompositionLocalAlignPair(char* seq1, char* seq2, int len1,
                            int len2, int alpha, int indel, int* submatrix, int* M);


/* added by Gary Benson 1/2003 */
int* GetCompositionSubstringMatchLengths(char* seq1, char* seq2, int len1, int len2,
																				 int limit, int ALPHABETSIZE);

/* added by Gary Benson 1/2003 */
void FreeCompositionAlignPair(COMPOSITIONALIGNPAIR* ap);

int*    CreateSubstitutionMatrix(int match,int mismatch);
char*   GetReverseComplement( char* original);


/************************************************************************
*                   - Macros, Constants, and Globals -
*************************************************************************/

/* contants, and macros */
#define ZERO            0
#define START           0
#define RIGHT           1
#define DOWN            2
#define DIAGONAL        3
#define WRAPRIGHT       4
#define WRAPDIAGONAL    5
#define LONGDIAGONAL    6
#define MINIMUM(a,b) (a<=b?a:b)
#define MAXIMUM(a,b) (a>=b?a:b)
#define min3switch(a,b,c) ((a<=b)?((a<=c)?1:3):((b<=c)?2:3))
#define max3switch(a,b,c) ((a>=b)?((a>=c)?1:3):((b>=c)?2:3))
#define max4switch(a,b,c,d) ((a>=b)?((a>=c)?((a>=d)?1:4):((c>=d)?3:4)):\
                            ((b>=c)?((b>=d)?2:4):((c>=d)?3:4)))


#define VLN             999999



/***********************************************************************
*                   - Function Definitions -
************************************************************************/


int*    CreateSubstitutionMatrix(int match,int mismatch)
{
    int *sm,*currint;
    int  i,j;

    /* allocate memory for matrix */
    sm = (int*) malloc(26*26*sizeof(int));
    if(sm==NULL) return NULL;

    /* generate 26x26 into matrix */
    for(i=1,currint=sm;i<=26;i++)
    {
        for(j=1;j<=26;j++,currint++)
        {
            if(i==j) *currint=match;
            else *currint=mismatch;
        }
    }

    /* return matrix */
    return sm;
}


/***************************************************************/
int* GetCompositionSubstringMatchLengths(char* seq1, char* seq2, int len1, int len2,
																				 int limit, int ALPHABETSIZE)
																				 /*
																				 Returns a two dimensional array (created and used as a one dimensional array) 
																				 M of size (len1+1)x(len2+1) which holds the length of the shortest 
																				 composition matching substrings (with zero < length <= limit) 
																				 ending at seq1[i] and seq2[j] or zero if no such matching substrings exist.
																				 Assumes the definition of a global ALPHABETSIZE.
																				 */
{
	
	__int64 MlenCheck; //to check for overflow
	int Mlen=0;
	COMPOSITIONVECTOR *Diffvectors, *dfi, *dfic;
  int *M, *Radsort1, *Radsort2, *Counts;
  int g,k,size,m,n,d,xstart,ystart,length;
	int minval,maxval,val,*rsi,*trs;
	int current,previous,match,i,j,substring_length;
  char *X, *Y;
  int *Index;	
	
	/* this check was put here because of a problem of overflow when calculating the Mlen
	and malloc() was just using the smaller number that came out as a result. 
	A very hard and pesky error. Gelfand. */
	MlenCheck=(__int64)(len1+1)*(__int64)(len2+1)*(__int64)(sizeof(int));
	if (MlenCheck>(__int64)INT_MAX) {
		return NULL;
	}
	
  /* allocate the substring composition match array */
	Mlen = (len1+1)*(len2+1);
	M = (int*) malloc(sizeof(int)*Mlen);
	if(M==NULL)
	{
		return NULL; /* in case memory allocation fails */
	}
	
  size=min(len1,len2);
  Diffvectors=(COMPOSITIONVECTOR *)malloc(sizeof(COMPOSITIONVECTOR)*(size+1));
	if(Diffvectors==NULL)
	{
		return NULL; /* in case memory allocation fails */
	}
  Radsort1=(int *)malloc(sizeof(int)*(size+1));
	if(Radsort1==NULL)
	{
		return NULL; /* in case memory allocation fails */
	}
  Radsort2=(int *)malloc(sizeof(int)*(size+1));
	if(Radsort2==NULL)
	{
		return NULL; /* in case memory allocation fails */
	}
  Counts=(int *)malloc(sizeof(int)*((2*size)+1));
	if(Counts==NULL)
	{
		return NULL; /* in case memory allocation fails */
	}
	Index=(int *)calloc(256,sizeof(int));
  if(Index==NULL){return NULL;}
  if(ALPHABETSIZE == 5)
  {
  Index['A']=0;
  Index['C']=1;
  Index['G']=2;
  Index['T']=3;
  Index['X']=4;
  }
  else
  {
  Index['A']=0;
  Index['B']=1;
  Index['C']=2;
  Index['D']=3;
  Index['E']=4;
  Index['F']=5;
  Index['G']=6;
  Index['H']=7;
  Index['I']=8;
  Index['J']=9;
  Index['K']=10;
  Index['L']=11;
  Index['M']=12;
  Index['N']=13;
  Index['O']=14;
  Index['P']=15;
  Index['Q']=16;
  Index['R']=17;
  Index['S']=18;
  Index['T']=19;
  Index['U']=20;
  Index['V']=21;
  Index['W']=22;
  Index['X']=23;
  Index['Y']=24;
  }
	/******/
	m=len1;
	n=len2;
	/* compute once for each diagonal of M */
	for(d=-m+1;d<=n-1;d++)
	{
		/* calculate starting locations in X and Y (strings start at 
		index zero) and the length of diagonal d */
		
		if(d<0)
		{
			xstart=0;
			ystart=-d;
			length=min(n,m+d);
		}
		else
		{
			xstart=d;
			ystart=0;
			length=min(n-d,m);
		}
    
		/* compute substring differences for diagonal d and store in Diffvectors */
		
		
		dfi=dfic=Diffvectors;
		Y=seq1+ystart;
		X=seq2+xstart;
		
		/* zeroth vector is all zeros */
		for(k=0;k<ALPHABETSIZE;k++)
			dfi->C[k]=0;
		
		/* step through the substrings computing the difference X-Y at each length */
		for(g=0;g<length;g++)
		{
			dfic++;
			/* current composition is equal to the previous composition */
			for(k=0;k<ALPHABETSIZE;k++)
			{
				dfic->C[k]=dfi->C[k];
			}
			/* and then adjusted for +X[g] and -Y[g] */
			dfic->C[Index[X[g]]]++;
			dfic->C[Index[Y[g]]]--;
			dfi++;
		}
		
	
		/* radixsort the composition vectors in Diffvectors */
		/* sort into Radsort1 using Radsort2 for swapping */
		/* length+1 is number of vectors to sort */
		/*********************/
		
		
		/* initialize Radsort1 to the order of the vectors in Diffvectors */
		for(g=0;g<=length;g++)
			Radsort1[g]=g;
		
		/* sort on each component starting with last */  
		for(k=ALPHABETSIZE-1;k>=0;k--)
		{
			
			/* get counts for counting sort */
			/* first zero counts */

			for(g=-length;g<=length;g++)
				Counts[g+length]=0;        /* +length so 0<= index <= 2*length */
			
			/* get counts in each component of vector */
			
			dfi=Diffvectors;
			minval=0;
			maxval=0;
			for(g=0;g<=length;g++)
			{
				val=dfi->C[k];
				Counts[val+length]++;
				if(val>maxval) maxval=val;
				if(val<minval) minval=val;
				dfi++;
			}
			
			/* make counts cumulative, use minval and maxval */
			
			for(g=minval+1;g<=maxval;g++)
				Counts[g+length]+=Counts[g+length-1];
			
				/* go through Diffvectors (last to first) and save new positions
			in Radsort2 */
			
			rsi=Radsort1+length;
			for(g=0;g<=length;g++)
			{
				val=Diffvectors[(*rsi)].C[k];
				Radsort2[Counts[val+length]-1]=(*rsi); /* subtract 1 because Radsort is zero based */
				Counts[val+length]--;
				rsi--;
			}
			
			/* swap Radsort2 into Radsort1 */
			trs=Radsort1;
			Radsort1=Radsort2;
			Radsort2=trs;
		}
		
		
		/**********************************/
		
		/* fill M */
		/* first fill zero row and zero column */
		for(g=0;g<=len2;g++)
			*(M+g)=0;
		for(g=1;g<=len1;g++)
			*(M+g*(len2+1))=0;
		
		/* find matching compositions in the sorted Diffvectors */
		/* length+1 is the number of vectors in the sorted list */
		
		rsi=Radsort1;
		previous=current=(*rsi);  /* previous not used the first time through the 
		loop */
		
		/* get location in M to save matching substring lengths */
		
		for(g=0;g<=length;g++)
		{
			/* first vector always has no previous one to look at so no match */
			if(g==0)
			{
				substring_length=0;
			}
			else
			{
				/* test if current and previous vectors match */
				/* but only if the matching substring would be no longer than limit */
				if((substring_length=current-previous)>limit) 
					substring_length=0; /* store a default length of zero */
				else
				{
					match=1;
					for(k=0;k<ALPHABETSIZE;k++)
					{
						if(Diffvectors[current].C[k]!=Diffvectors[previous].C[k])
						{
							match=0;
							break;
						}
					}
					
					/* if they match, compute matching substring length */
					if(match)
						substring_length=current-previous;
					else /* store a default length of zero */
						substring_length=0;
					/* previous=current; */ /* uncomment for longest substring match */
				}
				/* calculate substring ending indices */
			}
			if(d>=0)
			{
				j=current+d;
				i=current;
			}
			else
			{
				j=current;
				i=current-d;
			}
			
			/* store substring length at M[i][j] */
			*(M+i*(len2+1)+j)=substring_length;
			
			/* update pointers */
			previous=current;  /* uncomment for shortest substring */
			rsi++;
			current=(*rsi);
			/**********************************/
			/******/
		}
	}
	free(Index);
	free(Diffvectors);
	free(Radsort1);
	free(Radsort2);
	free(Counts);
	return(M);
}

/***********************************************************************/
COMPOSITIONALIGNPAIR* GetCompositionGlobalAlignPair(char* seq1, char* seq2, int len1,
                            int len2, int alpha, int indel, int* submatrix, int* M)
/*
Produces a composition alignment for two sequences.  The aligned pair has an 
extra string (matchboundaries) which marks the boundaries of the shortest composition matching substrings.
M[i][j] is a matrix (len1+1)x(len2+1) which holds the length of the shortest 
composition matching substrings (with zero< length < limit) ending at seq1[i] and seq2[j] 
or zero if no such matching substrings exist.
The limit is the longest of these short substring that can be included
in the alignment. 
M is computed with the function GetCompositionSubstringMatchLengths which includes 
the limit as a parameter.  This function must be passed M as a parameter.
In the alignment, a substring composition match is scored as the length of the 
substring times a matchscore which is assumed to be a constant.  Here it is set 
to the score of A vs A.
*/
{
    int     col, row,length;
    int     eb;             /* equals e+indel */
    int     fb;             /* equals f+indel */
    int     sa;             /* equals s+indel */
    int     stemp;          /* Gary Benson 10/13/03 */
    SD* S=NULL;
    int  Slen = 0;
    SD *Sp, *Srm1p, *Scm1p, *Srcm1p; /* temporary pointers */
	__int64 SlenCheck; //to check for overflow

    COMPOSITIONALIGNPAIR* ap;
    char *seq1pos,*seq2pos,*seq1sidepos,*seq2sidepos;

    int substringlen, matchscore, satype, i;
		char *mbpos;
		int *Mp;

		/* matchscore for substring composition match is the score of A vs A */
		matchscore=submatrix[26*('A'-'A')+('A'-'A')];

	/* this check was put here because of a problem of overflow when calculating the Slen
	and malloc() was just using the smaller number that came out as a result. A very hard and
	pesky error. Gelfand. */
	SlenCheck=(__int64)(len1+1)*(__int64)(len2+1)*(__int64)(sizeof(SD));
	if (SlenCheck>(__int64)INT_MAX) {
		return NULL;
	}

    /* allocate the alignment matrix */
    Slen = (len1+1)*(len2+1);
    S = (SD*) malloc(sizeof(SD)*Slen);
    if(S==NULL)
    {
        return NULL; /* in case memory allocation fails */
    }

	
    /* compute SDD Matrix */
    /* first element */
    Sp=S;
    Sp->score = 0;
    Sp->direction = START;

    /* across the top */
    for(col=1; col<=len2; col++)
    {
        Sp++;
        Sp->score = col*indel;
        Sp->direction = RIGHT;
    }

    /* Down the left side */
    Sp=S;
    for(row=1; row<=len1; row++)
    {
        Sp+=(len2+1);
        Sp->score = row*indel;
        Sp->direction = DOWN;
    }

    /* body of matrix */
    for(row=1;row<=len1;row++)
    {
        /* set all pointer around column one of current row */
        Sp      =   S+1+row*(len2+1);
        Scm1p   =   Sp-1;
        Srm1p   =   Scm1p-len2;
        Srcm1p  =   Srm1p-1;
				Mp			=		M+1+row*(len2+1);

        for(col=1;col<=len2;col++)
        {
            /* E */
            eb = Scm1p->score+indel;

            /* F */
            fb = Srm1p->score+indel;

            /* S */
						
						/* changed S computation, Gary Benson, 10/13/03 to allow comparison of mismatch
						and composition match greater than length 1 */
						
						/* calculate match or mismatch */
						
						sa = Srcm1p->score+submatrix[26*(seq1[row-1]-'A')+(seq2[col-1]-'A')];
						satype=DIAGONAL;
						
						/* calculate composition match and test against mismatch */
						if((substringlen=*Mp)>1)
						{
							stemp=(Sp-substringlen*(len2+2))->score+substringlen*matchscore;
                            if(stemp>sa) {
                                sa=stemp;
								satype=LONGDIAGONAL;
                            }
						}
						
		    switch(max3switch(sa,eb,fb))			
           {
                case 1:
                    Sp->score = sa;
                    Sp->direction = satype; /* DIAGONAL	or LONGDIAGONAL */
                    break;
                case 2:
                    Sp->score = eb;
                    Sp->direction = RIGHT;
                    break;
                case 3:
                    Sp->score = fb;
                    Sp->direction = DOWN;
                    break;
            }

            /* update pointers */
            Sp++;
            Srm1p++;
            Scm1p++;
            Srcm1p++;
						Mp++;
        }
    }

    /* get the length of the alignment */
    Sp = S+(len2+1)*len1+len2; /* set at last element */
		Mp = M+(len2+1)*len1+len2;

    /* find out how long the alignment is */
    length = 0;
    while(Sp->direction!=START)
    {
    /*     length ++; */
        switch(Sp->direction)
        {
            case START:
                 break;
            case DIAGONAL:
                 Sp = Sp -(len2+2);
								 Mp = Mp -(len2+2);
								 length++;
                 break;
						case LONGDIAGONAL:
							   substringlen=*Mp;
								 Sp=Sp-substringlen*(len2+2);
								 Mp=Mp-substringlen*(len2+2);
								 length+=substringlen;
								 break;
            case RIGHT:
                 Sp = Sp - 1;
								 Mp = Mp - 1;
								 length++;
                 break;
            case DOWN:
                 Sp = Sp - (len2+1);
								 Mp = Mp - (len2+1);
								 length++;
                 break;
        }
    }

    /*  allocate memory     */
    ap = (COMPOSITIONALIGNPAIR*) malloc(sizeof(COMPOSITIONALIGNPAIR));
    if(ap==NULL) return NULL;

    ap->sequence1side = (char*) malloc((length+1)*sizeof(char));
    if(ap->sequence1side==NULL) { free(ap); return NULL;}

    ap->sequence2side = (char*) malloc((length+1)*sizeof(char));
    if(ap->sequence2side==NULL)
    {
        free(ap->sequence1side);
        free(ap);
        return NULL;
    }
    ap->matchboundaries = (char*) malloc((length+1)*sizeof(char));
    if(ap->matchboundaries==NULL)
    {
			  free(ap->sequence2side);
        free(ap->sequence1side);
        free(ap);
        return NULL;
    }

    /***************************
    * trace back and fill in ap
    ****************************/
    ap->length = length;

    Sp = S+(len2+1)*len1+len2;      /* set at last element */
    Mp = M+(len2+1)*len1+len2;      /* set at last element */

    ap->score  = Sp->score;
    ap->sequence1end = len1;
    ap->sequence2end = len2;
    ap->sequence1start = 1;
    ap->sequence2start = 1;

    seq1pos = seq1+len1-1;    /* position last character */
    seq2pos = seq2+len2-1;
    seq1sidepos = ap->sequence1side+length;  /* position at termination */
    seq2sidepos = ap->sequence2side+length;
    *seq1sidepos=*seq2sidepos='\0'; /* terminate strings */
    seq1sidepos--;  /* move down to last char */
    seq2sidepos--;

    mbpos=ap->matchboundaries+length; /* matchboundaries string */
		*mbpos='\0';
		mbpos--;


    /* Sp has ben set to last element above */
    while(Sp->direction!=START)
    {
        switch(Sp->direction)
        {
            case DIAGONAL:
                 *seq1sidepos = *seq1pos;
                 *seq2sidepos = *seq2pos;
								 if(*seq1sidepos==*seq2sidepos) *mbpos='|';
								 else *mbpos=' ';
                 seq1pos--;seq2pos--;seq1sidepos--;seq2sidepos--;
								 mbpos--;
					       Sp = Sp -(len2+2);
                 Mp = Mp -(len2+2);
                 break;
            case LONGDIAGONAL:
							   substringlen=*Mp;
								 for(i=1;i<=substringlen;i++)
								 {
									 *seq1sidepos = *seq1pos;
									 *seq2sidepos = *seq2pos;
									 if(i==1) *mbpos='>';
									 else if (i==substringlen) *mbpos='<';
									 else *mbpos='-';
									 seq1pos--;seq2pos--;seq1sidepos--;seq2sidepos--;
									 mbpos--;
									 Sp = Sp -(len2+2);
                   Mp = Mp -(len2+2);
                 }
								 break;
            case RIGHT:
                 *seq1sidepos='-';
                 *seq2sidepos=*seq2pos;
                 *mbpos=' ';
								 seq1sidepos--;
                 seq2sidepos--;
                 seq2pos--;
								 mbpos--;
                 Sp = Sp - 1;
                 Mp = Mp - 1;
                 break;
            case DOWN:
                 *seq1sidepos=*seq1pos;
                 *seq2sidepos='-';
								 *mbpos=' ';
                 seq1sidepos--;
                 seq1pos--;
                 seq2sidepos--;
								 mbpos--;
                 Sp = Sp - (len2+1);
                 Mp = Mp - (len2+1);
                 break;
        }
    }

    /* free alignment matrix */
    free(S);
    

    return ap;
}

/***********************************************************************/
/***********************************************************************/
COMPOSITIONALIGNPAIR* GetCompositionLocalAlignPair(char* seq1, char* seq2, int len1,
                            int len2, int alpha, int indel, int* submatrix, int* M)
/*
Produces a local composition alignment for two sequences.  The aligned pair has an 
extra string (matchboundaries) which marks the boundaries of the shortest composition matching substrings.
M[i][j] is a matrix (len1+1)x(len2+1) which holds the length of the shortest 
composition matching substrings (with zero< length < limit) ending at seq1[i] and seq2[j] 
or zero if no such matching substrings exist.
The limit is the longest of these short substring that can be included
in the alignment. 
M is computed with the function GetCompositionSubstringMatchLengths which includes 
the limit as a parameter.  This function must be passed M as a parameter.
In the alignment, a substring composition match is scored as the length of the 
substring times a matchscore which is assumed to be a constant.  Here it is set 
to the score of A vs A.
*/
{
    int     col, row,length;
    int     eb;             /* equals e+indel */
    int     fb;             /* equals f+indel */
    int     sa;             /* equals s+indel */
    int			stemp;          /* Gary Benson 10/13/03 */
    SD* S=NULL;
    int  Slen = 0;
    SD *Sp, *Srm1p, *Scm1p, *Srcm1p; /* temporary pointers */
	__int64 SlenCheck; //to check for overflow

    COMPOSITIONALIGNPAIR* ap;
    char *seq1pos,*seq2pos,*seq1sidepos,*seq2sidepos;

    int substringlen, matchscore, satype, i;
		char *mbpos;
		int *Mp;
		int bestscore,bestscorerow,bestscorecol;


		/* matchscore for substring composition match is the score of A vs A */
		matchscore=submatrix[26*('A'-'A')+('A'-'A')];

	/* this check was put here because of a problem of overflow when calculating the Slen
	and malloc() was just using the smaller number that came out as a result. A very hard and
	pesky error. Gelfand. */
	SlenCheck=(__int64)(len1+1)*(__int64)(len2+1)*(__int64)(sizeof(SD));
	if (SlenCheck>(__int64)INT_MAX) {
		return NULL;
	}

    /* allocate the alignment matrix */
    Slen = (len1+1)*(len2+1);
    S = (SD*) malloc(sizeof(SD)*Slen);
    if(S==NULL)
    {
        return NULL; /* in case memory allocation fails */
    }

	
    /* compute SDD Matrix */
		/* set up for finding best score */
		bestscore=0;
		bestscorerow=0;
		bestscorecol=0;

    /* first element */
    Sp=S;
    Sp->score = 0;
    Sp->direction = START;

    /* across the top */
    for(col=1; col<=len2; col++)
    {
        Sp++;
        Sp->score = 0;
        Sp->direction = START;
    }

    /* Down the left side */
    Sp=S;
    for(row=1; row<=len1; row++)
    {
        Sp+=(len2+1);
        Sp->score = 0;
        Sp->direction = START;
    }

    /* body of matrix */
    for(row=1;row<=len1;row++)
    {
        /* set all pointer around column one of current row */
        Sp      =   S+1+row*(len2+1);
        Scm1p   =   Sp-1;
        Srm1p   =   Scm1p-len2;
        Srcm1p  =   Srm1p-1;
				Mp			=		M+1+row*(len2+1);

        for(col=1;col<=len2;col++)
        {
            /* E */
            eb = Scm1p->score+indel;

            /* F */
            fb = Srm1p->score+indel;

            /* S */
						
						/* changed S computation, Gary Benson, 10/13/03 to allow comparison of mismatch
						and composition match greater than length 1 */
						
						/* calculate match or mismatch */
						
						sa = Srcm1p->score+submatrix[26*(seq1[row-1]-'A')+(seq2[col-1]-'A')];
						satype=DIAGONAL;
						
						/* calculate composition match and test against mismatch */
						if((substringlen=*Mp)>1)
						{
							stemp=(Sp-substringlen*(len2+2))->score+substringlen*matchscore;
                            if(stemp>sa) {
                                sa=stemp;
								satype=LONGDIAGONAL;
                            }
						}
	
            switch(max4switch(0,sa,eb,fb))
            {
								case 1:
										Sp->score=0;
										Sp->direction=START;
										break;
                case 2:
                    Sp->score = sa;
                    Sp->direction = satype; /* DIAGONAL	or LONGDIAGONAL */
                    break;
                case 3:
                    Sp->score = eb;
                    Sp->direction = RIGHT;
                    break;
                case 4:
                    Sp->score = fb;
                    Sp->direction = DOWN;
                    break;
            }
						if(Sp->score>bestscore)
						{
							bestscore=Sp->score;
							bestscorerow=row;
							bestscorecol=col;
						}

            /* update pointers */
            Sp++;
            Srm1p++;
            Scm1p++;
            Srcm1p++;
						Mp++;
        }
    }

    /* get the length of the alignment */
    Sp = S+bestscorerow*(len2+1)+bestscorecol; /* set at bestscore element */
		Mp = M+bestscorerow*(len2+1)+bestscorecol;

    /* find out how long the alignment is */
    length = 0;
    while(Sp->direction!=START)
    {
        switch(Sp->direction)
        {
            case START:
                 break;
            case DIAGONAL:
                 Sp = Sp -(len2+2);
								 Mp = Mp -(len2+2);
								 length++;
                 break;
						case LONGDIAGONAL:
							   substringlen=*Mp;
								 Sp=Sp-substringlen*(len2+2);
								 Mp=Mp-substringlen*(len2+2);
								 length+=substringlen;
								 break;
            case RIGHT:
                 Sp = Sp - 1;
								 Mp = Mp - 1;
								 length++;
                 break;
            case DOWN:
                 Sp = Sp - (len2+1);
								 Mp = Mp - (len2+1);
								 length++;
                 break;
        }
    }

    /*  allocate memory     */
    ap = (COMPOSITIONALIGNPAIR*) malloc(sizeof(COMPOSITIONALIGNPAIR));
    if(ap==NULL) return NULL;

    ap->sequence1side = (char*) malloc((length+1)*sizeof(char));
    if(ap->sequence1side==NULL) { free(ap); return NULL;}

    ap->sequence2side = (char*) malloc((length+1)*sizeof(char));
    if(ap->sequence2side==NULL)
    {
        free(ap->sequence1side);
        free(ap);
        return NULL;
    }
    ap->matchboundaries = (char*) malloc((length+1)*sizeof(char));
    if(ap->matchboundaries==NULL)
    {
			  free(ap->sequence2side);
        free(ap->sequence1side);
        free(ap);
        return NULL;
    }

    /***************************
    * trace back and fill in ap
    ****************************/
    ap->length = length;

    Sp = S+bestscorerow*(len2+1)+bestscorecol;      /* set at bestscore element */
    Mp = M+bestscorerow*(len2+1)+bestscorecol;      /* set at bestscore element */

    ap->score  = Sp->score;
    ap->sequence1end = bestscorerow; /* was len1 */
    ap->sequence2end = bestscorecol; /* was len2 */

    seq1pos = seq1+bestscorerow-1;    /* position last character */
    seq2pos = seq2+bestscorecol-1;
    seq1sidepos = ap->sequence1side+length;  /* position at termination */
    seq2sidepos = ap->sequence2side+length;
    *seq1sidepos=*seq2sidepos='\0'; /* terminate strings */
    seq1sidepos--;  /* move down to last char */
    seq2sidepos--;

    mbpos=ap->matchboundaries+length; /* matchboundaries string */
		*mbpos='\0';
		mbpos--;


    /* Sp has been set to bestscore element above */
    while(Sp->direction!=START)
    {
        switch(Sp->direction)
        {
            case DIAGONAL:
                 *seq1sidepos = *seq1pos;
                 *seq2sidepos = *seq2pos;
								 if(*seq1sidepos==*seq2sidepos) *mbpos='|';
								 else *mbpos=' ';
                 seq1pos--;seq2pos--;seq1sidepos--;seq2sidepos--;
								 mbpos--;
					       Sp = Sp -(len2+2);
                 Mp = Mp -(len2+2);
                 break;
            case LONGDIAGONAL:
							   substringlen=*Mp;
								 for(i=1;i<=substringlen;i++)
								 {
									 *seq1sidepos = *seq1pos;
									 *seq2sidepos = *seq2pos;
									 if(i==1) *mbpos='>';
									 else if (i==substringlen) *mbpos='<';
									 else *mbpos='-';
									 seq1pos--;seq2pos--;seq1sidepos--;seq2sidepos--;
									 mbpos--;
									 Sp = Sp -(len2+2);
                   Mp = Mp -(len2+2);
                 }
								 break;
            case RIGHT:
                 *seq1sidepos='-';
                 *seq2sidepos=*seq2pos;
                 *mbpos=' ';
								 seq1sidepos--;
                 seq2sidepos--;
                 seq2pos--;
								 mbpos--;
                 Sp = Sp - 1;
                 Mp = Mp - 1;
                 break;
            case DOWN:
                 *seq1sidepos=*seq1pos;
                 *seq2sidepos='-';
								 *mbpos=' ';
                 seq1sidepos--;
                 seq1pos--;
                 seq2sidepos--;
								 mbpos--;
                 Sp = Sp - (len2+1);
                 Mp = Mp - (len2+1);
                 break;
        }
    }

		/* count characters in alignment to find start */
		length=0;
		for(row=0;row<ap->length;row++)
			if(ap->sequence1side[row]!='-') length++; 
    ap->sequence1start = ap->sequence1end-length+1;
		
		length=0;
		for(col=0;col<ap->length;col++)
			if(ap->sequence2side[col]!='-') length++; 
    ap->sequence2start = ap->sequence2end-length+1;
		
    /* free alignment matrix */
    free(S);
    

    return ap;
}

/***********************************************************************/


void FreeCompositionAlignPair(COMPOSITIONALIGNPAIR* ap)
{
    free(ap->sequence1side);
    free(ap->sequence2side);
		free(ap->matchboundaries);
    free(ap);
    return;
}

/********************************	init_complement_ascii				****************************/

char* init_complement_ascii(void)

{
  /* complement has 256 entries so that finding the entries for A, C, G and T which are complement ascii values */
  /* require no calculation */
  int i;
  char *Complementascii=(char *)calloc(256,sizeof(long int));
  
  if(Complementascii==NULL){ return NULL; }
  
  for (i=0; i<256; i++)
	Complementascii[i]=i;

  Complementascii['A']='T';
  Complementascii['C']='G';
  Complementascii['G']='C';
  Complementascii['T']='A';
  Complementascii['a']='t';
  Complementascii['c']='g';
  Complementascii['g']='c';
  Complementascii['t']='a';

  return Complementascii;
}

/***********************************************************************/
char*   GetReverseComplement(char* original)
{
    static char *Complementascii = NULL;            /* precomputed array used for the reverse complement function */
    char *sourceptr,*destinptr;
    int length;
    char *buffer;


    /* safety check */
	if (Complementascii==NULL) {
        Complementascii = init_complement_ascii();
        if ( NULL == Complementascii ) return NULL;
	}

    /* find out how long the string is */
    length = strlen(original);

    /* allocate the new string */
    buffer = (char*) malloc(sizeof(char)*(length+1));
    if(buffer==NULL)
    {
		return NULL;
    }

    /* reverse complement pattern */
    sourceptr = original;
    destinptr = buffer+length; /* position at termination */
    *destinptr = '\0'; /* terminate string */
    destinptr--;
    while(*sourceptr!='\0')
    {
		(*destinptr)=Complementascii[(*sourceptr)];
        destinptr--;
        sourceptr++;
    }

    /* return buffer containg reverse complement */
    return buffer;
}



   
#endif