/*-------------------------------------------------------------

Copyright (C) 2000 Peter Clote. 
All Rights Reserved.

Permission to use, copy, modify, and distribute this
software and its documentation for NON-COMMERCIAL purposes
and without fee is hereby granted provided that this
copyright notice appears in all copies.


THE AUTHOR AND PUBLISHER MAKE NO REPRESENTATIONS OR
WARRANTIES ABOUT THE SUITABILITY OF THE SOFTWARE, EITHER
EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE
IMPLIED WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
PARTICULAR PURPOSE, OR NON-INFRINGEMENT. THE AUTHORS
AND PUBLISHER SHALL NOT BE LIABLE FOR ANY DAMAGES SUFFERED
BY LICENSEE AS A RESULT OF USING, MODIFYING OR DISTRIBUTING
THIS SOFTWARE OR ITS DERIVATIVES.

-------------------------------------------------------------*/




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

	Program: multipleGotoh.c

 	Given input file of strings, all with same first
	and last letter, produces distance matrix. 
	This incorporates the correction to Nei, Waterman,
	but only computes the distances, rather than the
	alignment. The distances are INTEGER rather than
	float. In the program sequence
		phy0.c
		phy1.c
		phy2.c
	it is assumed that the distances are floats.


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

#include <stdio.h>
#include <stdlib.h>     /* prototype for calloc, malloc */
#include <string.h>     /* strcmp() prototype           */


#define N 10		/* upper bound on number strings */
#define MAXLEN 80       /* max length of string (name)  */
#define INF  30000
#define SPACE ' '
#define MIN2(x,y)   ( (x) <= (y) ? (x) : (y) )
#define MAX2(x,y)   ( (x) <= (y) ? (y) : (x) )
#define MIN(x,y,z)  MIN2(MIN2(x,y),z)
#define MISMATCH  3     /* cost of mismatched letters */
#define D(i,j)      d[(i)*(m+1) + (j)]
#define P(i,j)      p[(i)*(m+1) + (j)]
#define Q(i,j)      q[(i)*(m+1) + (j)]


#define U  2
#define V  2
#define W(i)    U*(i)+V       /* cost of  gap of length i, i>0 */
#define S(i)    s[(i)-1]
#define T(i)    t[(i)-1]
	/* D is integer matrix with n+1 rows and m+1 cols,
	   where n = length of s, m = length of t */
#define GAMMA(i,j) ( S(i) == T(j) ? 0 : MISMATCH )


typedef char String[MAXLEN];

void error ( char *s );     /* prototype of error function */

void error ( char *s )
{
	fprintf(stderr,"%s\n",s);
	exit(1);
}


/**************** MAIN BODY ********************/


main(int argc, char *argv[])
{
	String s,t;
	int distances[N][N];		/* holds distances between strings */
	int *d;
	int i,j,n,m;

	int *dist ( char *s, char *t );
	void display(int *, char *, char *);

	FILE *in;
	int num_strings;
	int num_lines(FILE*);
	String* strings, *initialize_array(FILE*, int);
	void display_array(String*, int);
	void display_distances( int d[N][N], int);
	void fancy_display_distances(int d[][], int );

	/*********  ERROR CHECKING in command line arguments *********/

	if (argc != 2)
		{
		fprintf(stderr,"Usage: %s filename.\n",argv[0]);
		exit(1);
     		}
	if ((in = fopen(argv[1],"r")) == NULL)
    		{ 
     		fprintf(stderr,"%s is not a readable file.\n", argv[1]);
     		exit(1);
     		}


	/*********  Compute size of array and initialize   ******/
   	num_strings =  num_lines(in);
   	strings = initialize_array(in,num_strings);
	if ( N < num_strings )
		{
     		fprintf(stderr,"Set constant N to %d\n",num_strings);
     		exit(1);
     		}

	for (i=0;i<num_strings;i++)
		{
		strcpy(s,strings[i]);
		n = strlen(s);
		for (j=i+1;j<num_strings;j++)
			{
			strcpy(t,strings[j]);
			m = strlen(t);
			d=dist(s,t);
			distances[i][j] = D(n,m);
			}
		}
	for (i=0;i<num_strings;i++)
		for (j=0;j<i;j++)
			distances[i][j] = distances[j][i];
	for (i=0;i<num_strings;i++)
		distances[i][i] = 0;
	display_distances(distances,num_strings);


}	/* end of main program */


int *dist ( char *s, char *t )
{
	int i,j,k,n = strlen(s), m = strlen(t);
	int *d,*p,*q;

	if ( s[0] != t[0] || s[n-1] != t[m-1] )
	  error("Strings don't have same beginning or ending character.");

	/******* allocate arrays *******************/
	d = (int *) calloc((n+1)*(m+1), sizeof(int));
	p = (int *) calloc((n+1)*(m+1), sizeof(int));
	q = (int *) calloc((n+1)*(m+1), sizeof(int));

	/******** initialize border values of matrices */
	D(0,0) = P(0,0) = Q(0,0) = 0;
	for (i=1;i<=n;i++)
		{
		D(i,0) =  W(i);
		Q(i,0) = INF;
		}
	for (j=1;j<=m;j++) 
		{
		D(0,j) = W(j);
		P(0,j) = INF;
		}

	/******** compute interior values using recurrence relation.
	 This loop is really of the form
				for k = 2 to n+m
					for i,j satisfying i+j = k
						.....
	*************************************************************/

	for (k=2;k<=n+m;k++)
		for (i=1;i<=MIN2(n,k-1);i++)
			{
			j = k-i;
			if ( j > m ) continue;     /* index out of bounds */
			P(i,j) = MIN2( D(i-1,j) + W(1), P(i-1,j) + U );
			Q(i,j) = MIN2( D(i,j-1) + W(1), Q(i,j-1) + U );
			D(i,j) = MIN( D(i-1,j-1) + GAMMA(i,j), P(i,j), Q(i,j) );
			}
	return(d);
}





void display ( int *d, char *s, char *t )
{
	/* display the distance matrix d(i,j) between strings
	   s,t where
			0<=i<=n = strlen(s)
	   and
			0<=j<=m = strlen(t) */

	int i,j, n = strlen(s), m = strlen(t);

	printf("Distance Matrix for alignment of %s with %s.\n\n",s,t);

	/****************** print letters of string t ***********/
	printf("%5c",SPACE);
	printf("%5c",SPACE);
	for (j=1;j<=m;j++) printf("%5c",T(j));
	printf("\n");


	/********************* print D(0,0),...,D(0,m) ************/
	printf("%5c",SPACE);
	for (j=0;j<=m;j++)
		printf("%5d",D(0,j));
	printf("\n");


	for (i=1;i<=n;i++)
		{
		printf("%5c",S(i));
		for (j=0;j<=m;j++)
			printf("%5d",D(i,j));
		printf("\n");
		}

}




void convert(char *p)
{
   /* fgets() keeps the \n marker, but when we compare
     strings using strcmp(), the presence of the marker
     may keep a match from being noticed. Hence we strip
     the \n marker */

   while (*p != '\n') ++p;
   *p = '\0';
}


String *initialize_array(FILE *in, int n)
{
   /* Assume that file contains strings, one per line */
   int i;
   String *strings;
   void convert(char *);

   strings = (String*) calloc(n,sizeof(String));
   for (i=0;i<n;++i)
      {
      fgets(strings[i],MAXLEN,in);
      convert(strings[i]);
      }
   fclose(in);
   return(strings);
}

int num_lines(FILE *in)
{
   int n=0,c;

   while ((c = fgetc(in)) != EOF)
     if (c == '\n') ++n;
   rewind(in);
   return(n);
}



/*---------------------------------
Procedure to output pairwise distances from
input file of strings. For debugging, use
the procedure
	fancy_display_distances
which prints row and column numbers.
-----------------------------------*/
void display_distances(int distances[N][N], int n)
{
	int i,j;

	for (i=0;i<n;i++)
		{
		for (j=0;j<n;j++)
			printf("%5d",distances[i][j]);
		printf("\n");
		}


}

/*---------------------------------
Procedure to output pairwise distances from
input file of strings. For debugging, use
the following procedure which prints row and column numbers.
-----------------------------------*/
void fancy_display_distances(int distances[N][N], int n)
{
	int i,j;

	
	/********************* print column numbers ****/
	printf("%8c",SPACE);
	for (j=0;j<n;j++)
		printf("%5d",j);
	printf("\n");
	printf("%8c",SPACE);
	for (j=0;j<n*5;j++) printf("-");
	printf("\n");



	for (i=0;i<n;i++)
		{
		printf("%5d | ",i);
		for (j=0;j<n;j++)
			printf("%5d",distances[i][j]);
		printf("\n");
		}


}