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

struct text {
	char *start;
	int length;
};

struct text slurp(char *filename) {
	FILE *file;
	struct text text;
	
	file = fopen(filename, "r");
	fseek(file, 0, SEEK_END);
	text.length = (int)ftell(file);
	fseek(file, 0, SEEK_SET);
	text.start = malloc(text.length);
	fread(text.start, 1, text.length, file);

	return text;
}

int *letter_frequency(struct text text) {
	char *p;
	char *e;
	int *freq;

	freq = calloc(128, sizeof(int));
	
	p = text.start;
	e = p + text.length;

	for (; p != e; ++p) {
		assert(*p >= 0);
		freq[*p]++;
	}

	return freq;
}

int *pair_frequency(struct text text) {
	char *p;
	char *e;
	int *freq;

	freq = calloc(128*128, sizeof(int));

	p = text.start;
	e = p + text.length - 1;

	for (; p != e; ++p) {
		int I = *p * 128 + p[1];
		freq[*p * 128 + p[1]]++;
	}

	return freq;
}

struct text index_from_letter_freq(int *freq) {
	struct text index;
	int i;
	char *p;

	index.length = 0;
	
	for (i=0; i<128; ++i)
		if (freq[i]) index.length++;

	index.start = malloc(index.length);

	for (i=0, p=index.start; i<128; ++i)
		if (freq[i]) *p++ = (char)i;

	return index;
}

static int *sort_index_freq;

int sort_index_cmp(const void *a, const void *b) {
	char A, B;
	A = *((char *)a);
	B = *((char *)b);
	if (sort_index_freq[A] > sort_index_freq[B]) return -1;
	if (sort_index_freq[A] == sort_index_freq[B]) return 0;
	return 1;
}

void sort_index(struct text index, int *freq) {
	sort_index_freq = freq;
	qsort(index.start, index.length, 1, sort_index_cmp);
}

char trnl(char c) {
	if (c == 10) return ' ';
	return c;
}

void print_letter_frequency(int *freq, struct text index) {
	int i;
	for (i=0; i<index.length; ++i) {
		int I = index.start[i];
		if (freq[I] == 0) continue;
		printf("%3d %c %8d\n", i, trnl((char)I), freq[I]);
	}
}

void print_pair_freq(int *freq, struct text index) {
	int i, j;
	int l;
	int max_freq;
	char *scale;
	int max_scale;

	l = index.length;

	max_freq = maximum(freq, l*l);
	
	scale = ".,,,-----oooooooo***************O";
	max_scale = strlen(scale);
	
	printf(" ");
	for (j=0; j<l; ++j) {
		printf("%c", trnl(index.start[j]));
	}
	
	printf("\n");

	for (i=0; i<l; ++i) {
		int I = index.start[i];
		printf("%c", trnl(I));
		for (j=0; j<index.length; ++j) {
			int J = index.start[j];
			int f = freq[l*i + j];
			char c = f == 0 ? ' ' : scale[(f * max_scale) / (max_freq + 1)];
			printf("%c", c);
		}
		printf("\n");
	}
}

int maximum(int *list, int length) {
	int max;
	int *end = list + length;
	
	max = 0;
	for (; list != end; ++list) {
		if (*list > max) max = *list;
	}

	return max;
}

struct before_after {
	double *before;
	double *after;
};

struct before_after before_after_probs(int *pair_freq, int n) {
	struct before_after probs;
	int i, j;

	probs.before = calloc(n*n, sizeof(double)); // given an i, prob that j follows
	probs.after = calloc(n*n, sizeof(double)); // given an i, prob that j was before
	
	for (i=0; i<n; ++i) {
		int total_before = 0;
		int total_after = 0;
		for (j=0; j<n; ++j) {
			total_before += pair_freq[i*n + j];
			total_after += pair_freq[j*n + i];
		}
		for (j=0; j<n; ++j) {
			probs.before[i*n + j] = 1.0 * pair_freq[i*n + j] / total_before;
			probs.after[i*n + j] = 1.0 * pair_freq[j*n + i] / total_after;
		}
	}

	return probs;
}

struct correlation {
	char i, j;
	double c;
};

int correlation_cmp(const void *a, const void *b) {
	const struct correlation *A = a;
	const struct correlation *B = b;
	if (A->c > B->c) return 1;
	if (A->c == B->c) return 0;
	return -1;
}

struct correlation *correlate(struct before_after probs, int n) {
	struct correlation *correlations;
	int i, j, k;
	correlations = calloc(n*n, sizeof(struct correlation));

	for (i=0; i<n; ++i) {
		for (j=0; j<n; ++j) {
			double total_dist = 0;
			for (k=0; k<n; ++k) {
				double d = probs.before[i*n+k] - probs.before[j*n+k];
				total_dist += d*d;
				d = probs.after[i*n+k] - probs.after[j*n+k];
				total_dist += d*d;
			}
			correlations[i*n + j].c = total_dist;
			correlations[i*n + j].i = (char)i;
			correlations[i*n + j].j = (char)j;
		}
	}

	return correlations;
}

int *compact_pairs(int *pair_freq, struct text index) {
	int i, j;
	int l = index.length;
	int *compact = calloc(l*l, sizeof(int));
	
	for (i=0; i<l; ++i) {
		int I = index.start[i];
		for (j=0; j<l; ++j) {
			int J = index.start[j];
			compact[i*l + j] = pair_freq[I*128 + J];
		}
	}

	return compact;
}

void print_correlations(struct correlation *c, int n, struct text index) {
	for (; n>0; --n, ++c) {
		char i = index.start[c->i];
		char j = index.start[c->j];
		printf("%c == %c : %7f\n", i, j, c->c);
	}
}

double maximum_correl(struct correlation *c, int n) {
	double max = 0;
	for (; n>0; --n, ++c) {
		if (c->c > max) max = c->c;
	}
	return max;
}

void print_correlations_table(struct correlation *correlations, struct text index) {
	int i, j;
	int l;
	double max_freq;
	char *scale;
	int max_scale;

	l = index.length;

	max_freq = maximum_correl(correlations, l*l);
	
	scale = " .,`/\\-=+x*XoO#%@";
	max_scale = strlen(scale);
	
	printf(" ");
	for (j=0; j<l; ++j) {
		printf("%c", trnl(index.start[j]));
	}
	
	printf("\n");

	for (i=0; i<l; ++i) {
		int I = index.start[i];
		printf("%c", trnl(I));
		for (j=0; j<index.length; ++j) {
			int J = index.start[j];
			double f = correlations[l*i + j].c;
			char c = f == 0 ? ' ' : scale[(int)((f * max_scale) / (max_freq + 0.0001))];
			printf("%c", c);
		}
		printf("\n");
	}
}

int main(int argc, char** argv) {
	struct text text;
	int *freq;
	struct text index;
	int *pair_freq;
	int n;
	struct before_after ba_probs;
	struct correlation *correlations;
	
	text = slurp("bible.txt");

	freq = letter_frequency(text);
	index = index_from_letter_freq(freq);
	n = index.length;
	index.start = " \n.!?:;,()-'0123456789AEIUYWVHFSZJCQKGTBPMNLORDaeiouywvhfszjcqkgtdbplrmnx";
//	sort_index(index, freq);
	pair_freq = pair_frequency(text);
	pair_freq = compact_pairs(pair_freq, index);
	ba_probs = before_after_probs(pair_freq, n);

	print_letter_frequency(freq, index);
	print_pair_freq(pair_freq, index);
	correlations = correlate(ba_probs, n);
	print_correlations_table(correlations, index);
	qsort(correlations, n*n, sizeof(struct correlation), correlation_cmp);
//	print_correlations(correlations, n*n, index);
}