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

#include <sys/types.h>
#include <sys/stat.h>

#define ARRAY 50
#define FILEBUF 1024

/* Macro's added by Patrick, mucho neato, aparently ifndef aren't POSIX
 * anymore */
#if !defined(MAX)
#define	MAX(A, B)	((A) > (B) ? (A) : (B))
#endif

#if !defined(MIN)
#define	MIN(A, B)	((A) < (B) ? (A) : (B))
#endif

/* Macro's added by Jon, allows checking of NULL pointers from memory
 * allocations without function calls, but how can you return values? */
#if !defined(MALLOC)
#define MALLOC(PTR, SIZE) do { PTR = malloc(SIZE); if (PTR == NULL) fatal_sys_error("malloc returned NULL"); } while (0)
#endif

#if !defined(REALLOC)
#define REALLOC(PTR, SIZE) do { PTR = realloc(PTR, SIZE); if (PTR == NULL && SIZE != 0) fatal_sys_error("realloc returned NULL"); } while (0)
#endif

/* Structure for sequences from FASTA files */
struct fasta {
	char *id;
	char *seq;
	int start;
	int end;
	int numres;
};

void fatal_error (char *message);
FILE * xfopen (const char *path, const char *mode);
void * xmalloc (size_t size);
void * xrealloc (void *ptr, size_t size);
int diff (char a, char b);
float pairwise (struct fasta *a, struct fasta *b);
void populate (struct fasta *a);
void check_length (struct fasta **array);
struct fasta ** read_fasta (FILE *fh);
char ** do_pairwise (struct fasta **array);
void free_fasta (struct fasta *);
int jsprintf (char *str, const char *const fmt, ...);

void
free_fasta (struct fasta *seq) {
	free(seq->id);
	free(seq->seq);
	free(seq);
}

void
fatal_error (char *message) {
	fprintf(stderr, "%s", message);
	exit(EXIT_FAILURE);
}

void
fatal_sys_error (char *message) {
	perror(message);
	exit(EXIT_FAILURE);
}

FILE *
xfopen (const char *path, const char *mode) {
	FILE *fh;
	fh = fopen(path, mode);
	if (fh == NULL)
		fatal_sys_error("fopen returned NULL");
	return fh;
}

void *
xmalloc (size_t size) {
	void *ptr;
	ptr = (void *) malloc(size);
	if (ptr == NULL)
		fatal_sys_error("malloc returned NULL");
	return ptr;
}

void *
xrealloc (void *ptr, size_t size) {
	ptr = (void *) realloc(ptr, size);
	if (ptr == NULL && size != 0)
		fatal_sys_error("realloc returned NULL");
	return ptr;
}

/* If the same resdidue return 1 else return 0 */
int 
diff (char a, char b) {
	if (a == b && a == '-')
		return 0;
	if (a == b)
		return 1;
	return 0;
}

/* Does the pairwise comparison
 * Compares the regions of the alignment that overlap, common gaps are
 * ignored, and the number of common residues divided by the length of
 * the shortest sequence is the identity.
 * */
float
pairwise (struct fasta *a, struct fasta *b) {
	float result;
	int start, end, numres, i, id = 0;
	
	/* If the sequences don't overlap then the seq ID is 0 */
	if (a->end < b->start || a->start > b->end)
		return 0;

	start = MIN( a->start, b->start );
	end = MAX( a->end, b->end );
	numres = MAX( a->numres, b->numres );

	for (i = start; i < end; i++)
		id += diff(a->seq[i], b->seq[i]);

	result = 100 * (float) id / (float) numres;

	return result;
}

/* "Populate" the rest of the fasta structure with information
 * Gets the start position of the sequence (first non-gap character)
 * Gets the end position (the last non-gap character)
 * Gets the number of residues in the sequence (that aren't gaps)
 * Gaps are represented by '-'
 * */
void
populate (struct fasta *a) {
	int i;
	int len = strlen(a->seq);

	a->numres = 0;
	
	for (i = 0; i <= len; i++) {
		if (a->seq[i] != '-') {
			a->start = i;
			break;
		}
	}
	for (i = len; i > 0; i--) {
		if (a->seq[i] != '-') {
			a->end = i;
			break;
		}
	}
	for (i = a->start; i < a->end; i++) {
		if (a->seq[i] != '-') {
			a->numres++;
		}
	}
}

/* Makes sure that all of the sequences are the same length */
void
check_length (struct fasta **array) {
	int i, length;

	if (array[0] != NULL)
		length = strlen(array[0]->seq);
	else {
		fatal_error("check_length() not passed an array of fasta structs\n");
	}

	for (i = 0; array[i] != NULL; i++) {
		if (length != (int) strlen(array[i]->seq)) {
			fatal_error("Not all of the sequences are the same length\n");
		}
	}
}

/* Reads in a FASTA file and returns an array of fasta structures */
struct fasta **
read_fasta (FILE *fh) {
	int i, j, k, c, filesize = 0;
	char *file;
	struct fasta **array;
	struct stat file_stat;

	array = (struct fasta **) malloc(ARRAY * sizeof(struct fasta *));

	/* Find the size of the file and allocate it */
	fstat(fileno(fh), &file_stat);
	if (! S_ISREG(file_stat.st_mode)) {
		fatal_error("Input file is not a regular file");
	}
	file = (char *) malloc(sizeof(char) * (file_stat.st_size + 1));

	/* Read in the file and null terminate it */
	while ((c = getc(fh)) != EOF) {
		file[filesize++] = c;
	}
	file[filesize] = '\0';

	/* Parse the FASTA file into an array of structures */
	for (i = 0, j = 0, k = 0; i < filesize; i++) {

		/* Begining of FASTA record */
		if (file[i] == '>') {
			if (j % ARRAY == 0)
				array = (struct fasta **) xrealloc(array, sizeof(struct fasta *) * (j + ARRAY));

			/* Initalise the structure */
			array[j] = (struct fasta *) xmalloc(sizeof(struct fasta));
			array[j]->id = (char *) xmalloc(sizeof(char));
			array[j]->seq = (char *) xmalloc(sizeof(char));
			array[j]->numres = 0;

			/* Read in the ID, terminated by a newline or NULL (end of string) */
			i++;
			while (file[i] != '\0' && file[i] != '\n') {
				if (k % FILEBUF == 0)
					array[j]->id = (char *) xrealloc(array[j]->id, sizeof(char) * (FILEBUF + k));
				array[j]->id[k++] = file[i++];
			}
			array[j]->id = (char *) xrealloc(array[j]->id, sizeof(char) * k + 1);
			array[j]->id[k] = '\0';
			k = 0;

			/* read in the sequence, terminated by a '>' */
			for (k = 0; file[i] != '\0' && file[i] != '>'; i++) {
				if (file[i] != '\n') {
					if (k % FILEBUF == 0)
						array[j]->seq = (char *) xrealloc(array[j]->seq, sizeof(char) * (FILEBUF + k));
					array[j]->seq[k++] = file[i];
				}
			}
			array[j]->seq = (char *) xrealloc(array[j]->seq, sizeof(char) * k + 1);
			array[j]->seq[k] = '\0';

			k = 0;
			i--;
			j++;
		}
	}

	free(file);
	array = (struct fasta **) xrealloc(array, sizeof(struct fasta *) * (j + 1));
	array[j] = NULL;

	/* find the start and end points for the alignments */
	for (i = 0; array[i] != NULL; i++)
		populate(array[i]); 
	return array;
}

char **
do_pairwise (struct fasta **array) {
	int j = 0, i = 0, offset = 0, size;
	char **output;

	/* Find the size of the output and preallocate the space in the output
	 * array */
	{
		while (array[i] != NULL)
			i++;

		/* Number of IDs, each ID, then the upper corner distances, then for
		 * NULL terminator */
		size = 1 + i + (i * i - i) / 2 + 1;

		output = (char **) xmalloc(sizeof(char *) * size);

		for (i = 0; i < size - 1; i++) {
			output[i] = (char *) xmalloc(sizeof(char) * FILEBUF);
			memset(output[i], 0x00, sizeof(char) * FILEBUF);
		}
		output[size - 1] = NULL;
	}

	/* The number of IDs */
	i = 0;
	while (array[i] != NULL)
		i++;
	jsprintf(output[offset++], "%i", i);

	/* If this statement is above the other one then it core dumps!? */

	/* Each ID */
	for (i = 0; array[i] != NULL; i++)
		jsprintf(output[offset++], "%s", array[i]->id);
		
	/* Upper corner pairwise comparisons */
	for (i = 0; array[i] != NULL; i++) {
		for (j = i + 1; array[j] != NULL; j++) {
			jsprintf(output[offset++], "%f", pairwise(array[i], array[j]));
		}
	}

	return output;
}

/* Grows buffer dynamically to make input fit inside it.
 * Requires malloc'd pointer to string.
 * Returns length of string.
 * */
int
jsprintf (char *str, const char *const fmt, ...) {
	int i, size;
	va_list ap;
	va_start(ap, fmt);

	str = (char *) realloc(str, sizeof(char) * FILEBUF);

	for (i = 1; (size = vsnprintf(str, sizeof(char) * FILEBUF * i, fmt, ap)) == -1;  i++) {
		if (i == 100)
			fatal_error("Reached 100 iterations of vsnprinft, it's probably blown up. Check input");
		str = (char *) xrealloc(str, sizeof(char) * FILEBUF * i);
	}

	str = (char *) xrealloc(str, sizeof(char) * size + 1);

	va_end(ap);
	return size;
}