/************************************************************************
*                       - function prototypes -
*************************************************************************/
/* added by Gary Benson 3/2003 */
COMPOSITIONALIGNPAIR* GetCompositionGlobalAlignPairSimpleMatchFunction(char* seq1, char* seq2, int len1,
                            int len2, int alpha, int indel, int* submatrix, int limit, 
														SUBSTRINGMATCHELEMENT* M);

/* added by Gary Benson 3/2003 */
COMPOSITIONALIGNPAIR* GetCompositionLocalAlignPairSimpleMatchFunction(char* seq1, char* seq2, int len1,
                            int len2, int alpha, int indel, int* submatrix, int limit, 
														SUBSTRINGMATCHELEMENT* M);

/* added by Gary Benson 8/2003 */
COMPOSITIONALIGNPAIR* GetCompositionPatternGlobalTextLocalAlignPairSimpleMatchFunction(char* seq1, char* seq2, int len1,
                            int len2, int alpha, int indel, int* submatrix, int limit, 
														SUBSTRINGMATCHELEMENT* M);
/*note pattern is sequence 2 */


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


/***********************************************************************/
COMPOSITIONALIGNPAIR* GetCompositionGlobalAlignPairSimpleMatchFunction
														(char* seq1, char* seq2, int len1,
                            int len2, int alpha, int indel, int* submatrix, int limit,
														SUBSTRINGMATCHELEMENT* 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 this version, the substring composition match has one of a number of simple scoring functions.
*/
{
    int     col, row,length;
    int     eb;             /* equals e+indel */
    int     fb;             /* equals f+indel */
    int     sa;             /* equals s+match or mismatch */
		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;
		SUBSTRINGMATCHELEMENT *Mp;


	/* check if the function has been set */
	if (pf==NULL) return NULL;


	/* 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;
						Mp->align_match_length=1;
						
						/* calculate composition match and test against mismatch */
						if((substringlen=Mp->substring_match_length)>1)
						{
							stemp=pf((Sp-substringlen*(len2+2))->score,row,col,
								substringlen,matchscore,matchscore);
							if(stemp>sa)
							{
                                sa=stemp;
								satype=LONGDIAGONAL;
								Mp->align_match_length=substringlen;
							}
						}
						
            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->align_match_length;
								 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->align_match_length;
								 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* GetCompositionLocalAlignPairSimpleMatchFunction(char* seq1, char* seq2, int len1,
                            int len2, int alpha, int indel, int* submatrix, int limit,
														SUBSTRINGMATCHELEMENT* 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 with a simple match function.  
*/
{
    int     col, row,length;
    int     eb;             /* equals e+indel */
    int     fb;             /* equals f+indel */
    int     sa;             /* equals s+match or mismatch */
		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;
		SUBSTRINGMATCHELEMENT *Mp;
		int bestscore,bestscorerow,bestscorecol;

	
	/* check if the function has been set */
	if (pf==NULL) return NULL;       
        

    /* 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;
				Mp->align_match_length=1;
				
				/* calculate composition match and test against mismatch */
				if((substringlen=Mp->substring_match_length)>1)
				{
					stemp=pf((Sp-substringlen*(len2+2))->score,row,col,
						substringlen,matchscore,matchscore);
					if(stemp>sa)
					{
                        sa=stemp;
						satype=LONGDIAGONAL;
						Mp->align_match_length=substringlen;
					}
				}
				
				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->align_match_length;
								 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->align_match_length;
								 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;
}


/***********************************************************************/
COMPOSITIONALIGNPAIR* GetCompositionPatternGlobalTextLocalAlignPairSimpleMatchFunction
														(char* seq1, char* seq2, int len1,
                            int len2, int alpha, int indel, int* submatrix, int limit,
														SUBSTRINGMATCHELEMENT* M)
/*
Produces a composition alignment for two sequences.  Global for the pattern, local for the text.
THE PATTERN MUST BE SEQUENCE 2. 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 this version, the substring composition match has one of a number of simple scoring functions.
*/
{
    int     col, row,length;
    int     eb;             /* equals e+indel */
    int     fb;             /* equals f+indel */
    int     sa;             /* equals s+match or mismatch */
		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;
		SUBSTRINGMATCHELEMENT *Mp;
    int bestscore,bestscorerow;


	/* check if the function has been set */
	if (pf==NULL) return NULL;


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


		/* set up to find bestscore in last column */
		bestscore=0;
		bestscorerow=0;

    /* 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 = 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;
						Mp->align_match_length=1;
						
						/* calculate composition match and test against mismatch */
						if((substringlen=Mp->substring_match_length)>1)
						{
							stemp=pf((Sp-substringlen*(len2+2))->score,row,col,
								substringlen,matchscore,matchscore);
							if(stemp>sa)
							{
                                sa=stemp;
								satype=LONGDIAGONAL;
								Mp->align_match_length=substringlen;
							}
						}
						
            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++;
        }

			  /* store bestscore in final column */
				if((Sp-1)->score>bestscore)
				{
					bestscore=(Sp-1)->score;
					bestscorerow=row;
				}
    }
    /* get the length of the alignment */
    Sp = S+(len2+1)*bestscorerow+len2; /* set at bestscore element in final column */
		Mp = M+(len2+1)*bestscorerow+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->align_match_length;
								 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)*bestscorerow+len2; /* set at bestscore element in final column */
		Mp = M+(len2+1)*bestscorerow+len2;

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

    seq1pos = seq1+bestscorerow-1;    /* position at bestscore character */
    seq2pos = seq2+len2-1;						/* position at end */
    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 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->align_match_length;
								 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 for sequence1 */
		length=0;
		for(row=0;row<ap->length;row++)
			if(ap->sequence1side[row]!='-') length++; 
    ap->sequence1start = ap->sequence1end-length+1;

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


    return ap;
}


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