#include <math.h>
#include <stdio.h>
#include <time.h>
#include <sys/time.h>
#include <cairo/cairo-pdf.h>
#include <X11/XKBlib.h>

#include "initcairo.h"
#include "triangle.h"
#include "linalg.h"

#define LOOP(i) for(int i = 0; i < 3; i++)

typedef struct {
	double x;
	double y;
} point_t;

void cartanMatrix(gsl_matrix *cartan, double a1, double a2, double a3, double s);
void initializeTriangleGenerators(gsl_matrix **gen, gsl_matrix *cartan);
point_t intersect(point_t a, point_t b, point_t c, point_t d);
void drawPoint(DrawingContext *ctx, point_t p);
void drawSegment(DrawingContext *ctx, point_t a, point_t b);
void drawLine(DrawingContext *ctx, point_t a, point_t b);
void drawImplicitLine(DrawingContext *ctx, double a, double b, double c);
void drawPolygon(DrawingContext *ctx, int sides, ...);
void drawAttractorConnection(DrawingContext *ctx, const char *word1, const char *word2);
void drawBox(DrawingContext *ctx, const char *word1, const char *word2);
void drawBoxStd(DrawingContext *ctx, const char *word, char base);
int compareAngle(const void *x, const void *y);
int computeLimitCurve();
void drawBoxes(DrawingContext *ctx);
void drawReflectors(DrawingContext *ctx);
void drawAttractors(DrawingContext *ctx);
void draw(DrawingContext *ctx);
void setup();
void destroy();
void print(DrawingContext *screen);
int processEvent(GraphicsInfo *info, XEvent *ev);

double parameter;
int n_group_elements;
double *limit_curve; // x, y, angle triples
int limit_curve_valid = 0;
groupelement_t *group;
workspace_t *ws;
gsl_matrix *gen[3];
gsl_matrix **matrices;
gsl_matrix *cartan, *cob, *coxeter, *coxeter_fixedpoints, *fixedpoints;

int show_boxes = 0;
int show_attractors = 0;
int show_reflectors = 0;
int show_limit = 1;
int limit_with_lines = 1;
int use_repelling = 1;

void cartanMatrix(gsl_matrix *cartan, double a1, double a2, double a3, double s)
{
	gsl_matrix_set(cartan, 0, 0, -2);
	gsl_matrix_set(cartan, 0, 1, 2*s*cos(a3));
	gsl_matrix_set(cartan, 0, 2, 2/s*cos(a2));

	gsl_matrix_set(cartan, 1, 0, 2/s*cos(a3));
	gsl_matrix_set(cartan, 1, 1, -2);
	gsl_matrix_set(cartan, 1, 2, 2*s*cos(a1));

	gsl_matrix_set(cartan, 2, 0, 2*s*cos(a2));
	gsl_matrix_set(cartan, 2, 1, 2/s*cos(a1));
	gsl_matrix_set(cartan, 2, 2, -2);
}

void initializeTriangleGenerators(gsl_matrix **gen, gsl_matrix *cartan)
{
	LOOP(i) gsl_matrix_set_identity(gen[i]);
	LOOP(i) LOOP(j) *gsl_matrix_ptr(gen[i], i, j) += gsl_matrix_get(cartan, i, j);
}

// intersect the lines a-b and c-d
point_t intersect(point_t a, point_t b, point_t c, point_t d)
{
	point_t res;
	double t = ((c.y-d.y)*(c.x-a.x) - (c.x-d.x)*(c.y-a.y)) / ((b.x-a.x)*(c.y-d.y) - (b.y-a.y)*(c.x-d.x));
	res.x = (1-t)*a.x + t*b.x;
	res.y = (1-t)*a.y + t*b.y;
	return res;
}

void drawPoint(DrawingContext *ctx, point_t p)
{
	cairo_t *C = ctx->cairo;

	cairo_save(C);
	cairo_move_to(C, p.x, p.y);
	cairo_close_path(C);
	cairo_set_line_cap(C, CAIRO_LINE_CAP_ROUND);
	cairo_set_line_width(C, 10.0/ctx->scalefactor);
	cairo_stroke(C);
	cairo_restore(C);
}

void drawSegment(DrawingContext *ctx, point_t a, point_t b)
{
	cairo_t *C = ctx->cairo;
	cairo_move_to(C, a.x, a.y);
	cairo_line_to(C, b.x, b.y);
	cairo_stroke(C);
}

void drawLine(DrawingContext *ctx, point_t a, point_t b)
{
	// not implemented
	point_t c, xplus, xminus;
	double normb;
	double scalbc;

	c.x = ctx->center_x;
	c.y = ctx->center_y;
	b.x -= a.x;
	b.y -= a.y;
	c.x -= a.x;
	c.y -= a.y;
	normb = sqrt(b.x*b.x + b.y*b.y);
	b.x /= normb;
	b.y /= normb;
	scalbc = b.x*c.x + b.y*c.y;
	xplus.x = a.x + (scalbc + ctx->radius)*b.x;
	xplus.y = a.y + (scalbc + ctx->radius)*b.y;
	xminus.x = a.x + (scalbc - ctx->radius)*b.x;
	xminus.y = a.y + (scalbc - ctx->radius)*b.y;

	drawSegment(ctx, xminus, xplus);
}

void drawImplicitLine(DrawingContext *ctx, double a, double b, double c)
{
	// not implemented
	double norm, lambda;
	point_t m, s, xminus, xplus;
	m.x = ctx->center_x;
	m.y = ctx->center_y;
	lambda = (a*m.x + b*m.y + c)/(a*a + b*b);
	s.x = m.x - lambda*a;
	s.y = m.y - lambda*b;
	norm = sqrt(a*a + b*b);
	xminus.x = s.x - ctx->radius * b / norm;
	xminus.y = s.y + ctx->radius * a / norm;
	xplus.x = s.x + ctx->radius * b / norm;
	xplus.y = s.y - ctx->radius * a / norm;

	drawSegment(ctx, xminus, xplus);
}

void drawPolygon(DrawingContext *ctx, int sides, ...)
{
	va_list args;
	point_t first, prev, current;

	va_start(args, sides);

	first = va_arg(args, point_t);
	current = first;
	for(int i = 0; i < sides-1; i++) {
		prev = current;
		current = va_arg(args, point_t);
		drawSegment(ctx, prev, current);
	}
	drawSegment(ctx, current, first);

	va_end(args);
}

int fixedPoints(const char *word, point_t *out)
{
	gsl_matrix *tmp = ws->stack[10];
	gsl_matrix *ev = ws->stack[11];

	gsl_matrix_set_identity(tmp);
	for(int i = 0; i < strlen(word); i++) {
		if(word[i] == ' ')
			continue;
		multiply_right(tmp, gen[word[i]-'a'], ws);
	}
	int count = real_eigenvectors(tmp, ev, ws);
	multiply_left(cob, ev, ws);

	LOOP(i) out[i].x = gsl_matrix_get(ev, 0, i) / gsl_matrix_get(ev, 2, i);
	LOOP(i) out[i].y = gsl_matrix_get(ev, 1, i) / gsl_matrix_get(ev, 2, i);

	return count;
}

void drawAttractorConnection(DrawingContext *ctx, const char *word1, const char *word2)
{
	point_t p1[3], p2[3];

	fixedPoints(word1, p1);
	fixedPoints(word2, p2);

	drawSegment(ctx, p1[0], p2[0]);
}

void drawTriangle(DrawingContext *ctx, const char *word)
{
	point_t p[3];

	fixedPoints(word, p);
	drawPolygon(ctx, 3, p[0], p[1], p[2]);
}

void drawBox(DrawingContext *ctx, const char *word1, const char *word2)
{
	point_t p[2][3],i[2];

	fixedPoints(word1, p[0]);
	fixedPoints(word2, p[1]);

	// intersect attracting line with neutral line of the other element
	for(int j = 0; j < 2; j++)
		i[j] = intersect(p[j%2][0],p[j%2][1],p[(j+1)%2][0],p[(j+1)%2][2]);

	drawPolygon(ctx, 4, p[0][0], i[0], p[1][0], i[1]);
}

void drawBoxStd(DrawingContext *ctx, const char *word, char base)
{
	char word1[100];
	char word2[100];

	int len = strlen(word);
	if(len*2 + 4 > 100)
		return;

	for(int i = 0; i < len; i++) {
		word1[i] = word1[2*len+2-i] = word[i];
		word2[i] = word2[2*len+2-i] = word[i];
	}
	word1[2*len+3] = 0;
	word2[2*len+3] = 0;

	LOOP(i) word1[len+i] = (base-'A'+6+i+1)%3+'a';
	LOOP(i) word2[len+i] = (base-'A'+6-i-1)%3+'a';

//	printf("Words: %s %s\n", word1, word2);

	drawBox(ctx, word1, word2);
}

int compareAngle(const void *x, const void *y)
{
	return ((double*)x)[2] > ((double*)y)[2] ? 1 : -1;
}

int computeLimitCurve()
{
	int p = 5, q = 5, r = 5;
	int k = 2, l = 2, m = 2;

	int column = use_repelling ? 2 : 0;

	generate_triangle_group(group, n_group_elements, p, q, r);

	// do first in the Fuchsian case to get the angles
	cartanMatrix(cartan, M_PI/p, M_PI/q, M_PI/r, 1.0);
	initializeTriangleGenerators(gen, cartan);
	gsl_matrix_set_identity(matrices[0]);
	for(int i = 1; i < n_group_elements; i++)
		multiply(matrices[group[i].parent->id], gen[group[i].letter], matrices[i]);
	diagonalize_symmetric_form(cartan, cob, ws);
	multiply_many(ws, coxeter, 3, gen[0], gen[1], gen[2]);
	int ev_count_fuchsian = real_eigenvectors(coxeter, coxeter_fixedpoints, ws);

	if(ev_count_fuchsian != 3)
		return 0;

	for(int i = 0; i < n_group_elements; i++) {
		multiply_many(ws, fixedpoints, 3, cob, matrices[i], coxeter_fixedpoints);
		limit_curve[3*i+2] = atan2(
			gsl_matrix_get(fixedpoints, 0, column)/gsl_matrix_get(fixedpoints, 2, column),
			gsl_matrix_get(fixedpoints, 1, column)/gsl_matrix_get(fixedpoints, 2, column));
	}

	// now to it again to calculate x and y coordinates
	cartanMatrix(cartan, M_PI*k/p, M_PI*l/q, M_PI*m/r, parameter);
	initializeTriangleGenerators(gen, cartan);
	gsl_matrix_set_identity(matrices[0]);
	for(int i = 1; i < n_group_elements; i++)
		multiply(matrices[group[i].parent->id], gen[group[i].letter], matrices[i]);
	gsl_matrix_memcpy(cob, cartan); // is this a good choice of basis
//	gsl_matrix_set_identity(cob);
	multiply_many(ws, coxeter, 3, gen[0], gen[1], gen[2]);
	int ev_count = real_eigenvectors(coxeter, coxeter_fixedpoints, ws);

	if(ev_count == 1)
		column = 0;
	if(ev_count == 0)
		return 0;

	for(int i = 0; i < n_group_elements; i++) {
		multiply_many(ws, fixedpoints, 3, cob, matrices[i], coxeter_fixedpoints);

		limit_curve[3*i] = gsl_matrix_get(fixedpoints, 0, column)/gsl_matrix_get(fixedpoints, 2, column);
		limit_curve[3*i+1] = gsl_matrix_get(fixedpoints, 1, column)/gsl_matrix_get(fixedpoints, 2, column);
	}

	qsort(limit_curve, n_group_elements, 3*sizeof(double), compareAngle);

	return 1;
}

void drawBoxes(DrawingContext *ctx)
{
	cairo_t *C = ctx->cairo;

	cairo_set_source_rgb(C, 1, 0, 0);
	drawTriangle(ctx, "abc");
	cairo_set_source_rgb(C, 0, 0, 1);
	drawTriangle(ctx, "aca abc aca");
	drawTriangle(ctx, "acac abc caca");
	drawTriangle(ctx, "acaca abc acaca");
	cairo_set_source_rgb(C, 0, 0.8, 0);
	drawTriangle(ctx, "cac abc cac");
	drawTriangle(ctx, "caca abc acac");
	drawTriangle(ctx, "cacac abc cacac");

	/*
	cairo_set_source_rgb(C, 1, 0, 0);
	drawBoxStd(ctx, "c", 'C');
	drawBoxStd(ctx, "", 'B');
	drawBoxStd(ctx, "a", 'A');
	drawBoxStd(ctx, "", 'C');
	drawBoxStd(ctx, "b", 'B');

	cairo_set_source_rgb(C, 0, 0, 0);
	drawBoxStd(ctx, "ca", 'A');
	drawBoxStd(ctx, "cac", 'C');
	drawBoxStd(ctx, "caca", 'A');
	drawBoxStd(ctx, "acac", 'C');
	drawBoxStd(ctx, "aca", 'A');
	drawBoxStd(ctx, "ac", 'C');

	drawBoxStd(ctx, "aca cb", 'B');
	drawBoxStd(ctx, "aca cbc", 'C');
	drawBoxStd(ctx, "aca cbcb", 'B');
	drawBoxStd(ctx, "aca bcbc", 'C');
	drawBoxStd(ctx, "aca bcb", 'B');
	drawBoxStd(ctx, "aca bc", 'C');

	drawBoxStd(ctx, "caca cb", 'B');
	drawBoxStd(ctx, "caca cbc", 'C');
	drawBoxStd(ctx, "caca cbcb", 'B');
	drawBoxStd(ctx, "caca bcbc", 'C');
	drawBoxStd(ctx, "caca bcb", 'B');
	drawBoxStd(ctx, "caca bc", 'C');

	cairo_set_source_rgb(C, 1, 0, 1);
	drawBoxStd(ctx, "ca bc", 'C');
	drawBoxStd(ctx, "ca bcb", 'B');
	drawBoxStd(ctx, "ca bcbc", 'C');
	drawBoxStd(ctx, "ca cbcb", 'B');
	drawBoxStd(ctx, "ca cbc", 'C');
	drawBoxStd(ctx, "ca cb", 'B');

	cairo_set_source_rgb(C, 0, 1, 0);
//	drawBoxStd(ctx, "ca bc", 'C');
	drawBoxStd(ctx, "cabc ba", 'A');
	drawBoxStd(ctx, "cabc bab", 'B');
	drawBoxStd(ctx, "cabc baba", 'A');
	drawBoxStd(ctx, "cabc abab", 'B');
	drawBoxStd(ctx, "cabc aba", 'A');
	drawBoxStd(ctx, "cabc ab", 'B');
	*/
}

void drawReflectors(DrawingContext *ctx)
{
	gsl_matrix *tmp = ws->stack[10];
	gsl_matrix *tmp2 = ws->stack[11];
	point_t p[3];

	// points
	cairo_set_source_rgb(ctx->cairo, 0, 0, 0);
	LOOP(i) p[i].x = gsl_matrix_get(cob, 0, i) / gsl_matrix_get(cob, 2, i);
	LOOP(i) p[i].y = gsl_matrix_get(cob, 1, i) / gsl_matrix_get(cob, 2, i);
	LOOP(i) drawPoint(ctx, p[i]);

	// lines
	invert(cob, tmp2, ws);
	multiply(cartan, tmp2, tmp);
	LOOP(i) drawImplicitLine(ctx, gsl_matrix_get(tmp, i, 0), gsl_matrix_get(tmp, i, 1), gsl_matrix_get(tmp, i, 2));
}

void drawAttractors(DrawingContext *ctx)
{
	point_t p[3][3];

	fixedPoints("abc", p[0]);
	fixedPoints("bca", p[1]);
	fixedPoints("cab", p[2]);

	double color[3][3] = {{1,0,0},{0,0.7,0},{0,0,1}};

	LOOP(i) LOOP(j) {
		cairo_set_source_rgb(ctx->cairo, color[i][0], color[i][1], color[i][2]);
		drawPoint(ctx, p[i][j]);
	}

	LOOP(i) LOOP(j) {
		cairo_set_source_rgb(ctx->cairo, color[i][0], color[i][1], color[i][2]);
		drawLine(ctx, p[i][j], p[i][(j+1)%3]);
	}
}

void draw(DrawingContext *ctx)
{
	cairo_t *C = ctx->cairo;

	cairo_set_source_rgb(C, 1, 1, 1);
	cairo_paint(C);

	cairo_set_matrix(C, &ctx->matrix);

	// defaults; use save/restore whenever these are changed
	cairo_set_line_width(C, 1.0/ctx->scalefactor);
	cairo_set_font_size(C, 16);
	cairo_set_line_join(C, CAIRO_LINE_JOIN_BEVEL);
	cairo_set_line_cap(C, CAIRO_LINE_CAP_ROUND);

	if(limit_curve_valid) {
		if(show_limit) {
			cairo_save(C);

			int previous_inside = 0;
			for(int i = 0; i < n_group_elements; i++) {
				double x = limit_curve[3*i];
				double y = limit_curve[3*i+1];

				cairo_user_to_device(C, &x, &y);

				if(-x < ctx->width && x < 3*ctx->width && -y < ctx->height && y < 3*ctx->height) {
					if(limit_with_lines) {
						if(!previous_inside)
							cairo_move_to(C, limit_curve[3*i], limit_curve[3*i+1]);
						else
							cairo_line_to(C, limit_curve[3*i], limit_curve[3*i+1]);
					} else {
						cairo_move_to(C, limit_curve[3*i], limit_curve[3*i+1]);
						cairo_close_path(C);
					}
					previous_inside = 1;
				} else {
					previous_inside = 0;
				}
			}

			if(!limit_with_lines) { // draw dots instead of lines
				cairo_set_line_cap(C, CAIRO_LINE_CAP_ROUND);
				cairo_set_line_width(C, 3.0/ctx->scalefactor);
			}

			cairo_set_source_rgb(C, 0, 0, 0);
			cairo_stroke(C);

			cairo_restore(C);
		}

		if(show_boxes)
			drawBoxes(ctx);

		if(show_attractors)
			drawAttractors(ctx);

		if(show_reflectors)
			drawReflectors(ctx);
	}

	cairo_identity_matrix(C);

	cairo_move_to(C, 15, 30);
	cairo_set_source_rgb(C, 0, 0, 0);
	char buf[100];
	sprintf(buf, "t = exp(%.8f) = %.8f", log(parameter), parameter);
	cairo_show_text(C, buf);

	cairo_surface_flush(cairo_get_target(C));
}

void setup()
{
	n_group_elements = 50000;
	limit_curve = malloc(3*n_group_elements*sizeof(double));
	group = malloc(n_group_elements*sizeof(groupelement_t));
	ws = workspace_alloc(3);

	LOOP(i) gen[i] = gsl_matrix_alloc(3, 3);
	matrices = malloc(n_group_elements*sizeof(gsl_matrix*));
	for(int i = 0; i < n_group_elements; i++)
		matrices[i] = gsl_matrix_alloc(3, 3);
	cartan = gsl_matrix_alloc(3, 3);
	cob = gsl_matrix_alloc(3, 3);
	coxeter = gsl_matrix_alloc(3, 3);
	coxeter_fixedpoints = gsl_matrix_alloc(3, 3);
	fixedpoints = gsl_matrix_alloc(3, 3);
}

void destroy()
{
	free(limit_curve);
	free(group);
	workspace_free(ws);

	LOOP(i) gsl_matrix_free(gen[i]);
	for(int i = 0; i < n_group_elements; i++)
		gsl_matrix_free(matrices[i]);
	free(matrices);
	gsl_matrix_free(cartan);
	gsl_matrix_free(cob);
	gsl_matrix_free(coxeter);
	gsl_matrix_free(coxeter_fixedpoints);
	gsl_matrix_free(fixedpoints);
}

void print(DrawingContext *screen)
{
	DrawingContext file;
	cairo_surface_t *surface;

	char filename[100];
	time_t t = time(NULL);
	strftime(filename, sizeof(filename), "screenshot_%Y%m%d_%H%M%S.pdf", localtime(&t));

	file.width = screen->width;
	file.height = screen->width / sqrt(2.0);
	file.matrix = screen->matrix;
	file.matrix.y0 += ((double)file.height - (double)screen->height) / 2.0; // recenter vertically
	updateDimensions(&file);

	surface = cairo_pdf_surface_create(filename, (double)file.width, (double)file.height);
    file.cairo = cairo_create(surface);

	draw(&file);

	cairo_destroy(file.cairo);
	cairo_surface_destroy(surface);

	printf("Wrote sceenshot to file: %s\n", filename);
}

int processEvent(GraphicsInfo *info, XEvent *ev)
{
	int state;
	unsigned long key;

	switch(ev->type) {

	case KeyPress:
		state = ev->xkey.state & (ShiftMask | LockMask | ControlMask);
		key = XkbKeycodeToKeysym(ev->xkey.display, ev->xkey.keycode, 0, !!(state & ShiftMask));
		printf("Key pressed: %ld\n", key);

		switch(key) {
		case XK_Down:
			parameter /= exp(0.002);
			limit_curve_valid = computeLimitCurve();
			break;
		case XK_Up:
			parameter *= exp(0.002);
			limit_curve_valid = computeLimitCurve();
			break;
		case XK_Left:
			parameter /= exp(0.0001);
			limit_curve_valid = computeLimitCurve();
			break;
		case XK_Right:
			parameter *= exp(0.0001);
			limit_curve_valid = computeLimitCurve();
			break;
		case XK_Page_Down:
			parameter /= exp(0.2);
			limit_curve_valid = computeLimitCurve();
			break;
		case XK_Page_Up:
			parameter *= exp(0.2);
			limit_curve_valid = computeLimitCurve();
			break;
		case ' ':
			parameter = 2.890053638;
			limit_curve_valid = computeLimitCurve();
			break;
		case XK_Return:
			parameter = 2.76375163;
			limit_curve_valid = computeLimitCurve();
			break;
		case 'm':
			printf("info->context->matrix.xx = %f;\n", info->context->matrix.xx);
			printf("info->context->matrix.xy = %f;\n", info->context->matrix.xy);
			printf("info->context->matrix.x0 = %f;\n", info->context->matrix.x0);
			printf("info->context->matrix.yx = %f;\n", info->context->matrix.yx);
			printf("info->context->matrix.yy = %f;\n", info->context->matrix.yy);
			printf("info->context->matrix.y0 = %f;\n", info->context->matrix.y0);
			break;
		case 'b':
			show_boxes = !show_boxes;
			break;
		case 'a':
			show_attractors = !show_attractors;
			break;
		case 'r':
			show_reflectors = !show_reflectors;
			break;
		case 'L':
			limit_with_lines = !limit_with_lines;
			break;
		case 'l':
			show_limit = !show_limit;
			break;
		case 'p':
			print(info->context);
			break;
		case 'f':
			use_repelling = !use_repelling;
			limit_curve_valid = computeLimitCurve();
		}

		return STATUS_REDRAW;
	}

	return STATUS_NOTHING;
}

int main()
{
	GraphicsInfo *info;

	parameter = 3.0;

	setup();
	limit_curve_valid = computeLimitCurve();

	info = initCairo(0, KeyPressMask, 200, 200, "Triangle group");
	if(!info)
		return 1;

	DrawingContext *ctx = info->context;

	ctx->matrix.xx = 837.930824;
	ctx->matrix.xy = -712.651341;
	ctx->matrix.x0 = 180.427716;
	ctx->matrix.yx = 712.651341;
	ctx->matrix.yy = 837.930824;
	ctx->matrix.y0 = 1412.553240;

	/*
	info->context->matrix.xx = 1636.583641;
	info->context->matrix.xy = -1391.897150;
	info->context->matrix.x0 = 975.772883;
	info->context->matrix.yx = 1391.897150;
	info->context->matrix.yy = 1636.583641;
	info->context->matrix.y0 = 2446.174297;
	*/

	updateDimensions(ctx);

	startTimer(info);

	while(1) {
		int result = checkEvents(info, processEvent, NULL);
		if(result == STATUS_QUIT)
			return 0;
		else if(result == STATUS_REDRAW) {
			struct timeval current_time;
			double start_time, intermediate_time, end_time;
			gettimeofday(&current_time, 0);
			start_time = current_time.tv_sec + current_time.tv_usec*1e-6;

			updateDimensions(info->context);
			draw(info->context);

			gettimeofday(&current_time, 0);
			intermediate_time = current_time.tv_sec + current_time.tv_usec*1e-6;

			cairo_set_source_surface(info->front_context, info->buffer_surface, 0, 0);
			cairo_paint(info->front_context);

			gettimeofday(&current_time, 0);
			end_time = current_time.tv_sec + current_time.tv_usec*1e-6;
			printf("drawing finished in %.2f milliseconds, of which %.2f milliseconds were buffer switching\n", (end_time - start_time) * 1000, (end_time - intermediate_time) * 1000);
		}
		waitUpdateTimer(info);
	}

	destroyCairo(info);

	destroy();

	return 0;
}