/*	a2serial.c

	2011 by Hank Dietz for Fall EE599/EE699
*/

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

/*	Maximum number of points */
#define	MAXLEN	(16*1024)

float	weight[MAXLEN][MAXLEN];
int	point[MAXLEN];
int	color[MAXLEN];
int	xdim = 0;
int	ydim = 0;
int	genelen = 0;
int	colors = 2;

#define	MAXMAP	(1024*1024)
unsigned char pixels[MAXMAP];

#define	MAXPOP	1024

typedef struct {
	float	prio[MAXLEN];
} geneome;

typedef	struct {
	geneome	g;
	float	v;
} individual;

individual	p[MAXPOP];


static int	nextc = ' ';

inline static int
advance(void)
{
	if (nextc != EOF) nextc = getchar();
	return(nextc);
}

inline static void
eat_space(void)
{
	while ((nextc == ' ') ||
	       (nextc == '\t') ||
	       (nextc == '\n') ||
	       (nextc == '\r') ||
	       (nextc == '#')) {
		if (nextc == '#') {
			while ((nextc != '\n') &&
			       (nextc != EOF)) {
				advance();
			}
		}
		advance();
	}
}

inline static int
get_number(void)
{
	register int n = 0;

	while ((nextc >= '0') && (nextc <= '9')) {
		n = (n * 10) + nextc - '0';
		advance();
	}

	return(n);
}

extern	float sqrtf(float);

inline static inline float
dist(register int a,
register int b)
{
	/* Distance: lower is better, nearer */
	register float dx = (a % xdim) - (b % xdim);
	register float dy = ((a / xdim) % ydim) - ((b / xdim) % ydim);
	dx *= dx;
	dy *= dy;
	return(sqrtf(dx + dy));
}

read_P5(void)
{
	/* Read and interpret a P5 file */
	register int i, j, xy, maxval;

	/* Read the header */
	if (advance() != 'P') {
badP5:
		fprintf(stderr,
			"Error: Input is not a usable P5 file\n");
		exit(1);
	}
	if (advance() != '5') goto badP5;
	advance();
	eat_space();
	xdim = get_number();
	if (xdim < 2) goto badP5;
	eat_space();
	ydim = get_number();
	if (ydim < 2) goto badP5;
	if ((xdim * ydim) > MAXMAP) goto badP5;
	eat_space();
	maxval = get_number();
	if ((maxval < 1) || (maxval > 255)) goto badP5;

	/* Read the pixel map counting points & saving offsets */
	genelen = 0;
	xy = xdim * ydim;
	for (i=0; i<xy; ++i) {
		if ((pixels[i] = advance()) != 0) {
			point[genelen] = i;
			if (++genelen >= MAXLEN) goto badP5;
		}
	}

	/* Build the weight table */
	for (i=0; i<genelen; ++i) {
		for (j=0; j<genelen; ++j) {
			weight[j][i] =
			weight[i][j] = dist(point[i], point[j]);
		}
	}
}

typedef	struct {
	int	index;
	float	key;
} eval_t;

static	eval_t	evals[MAXLEN];

float
eval(register geneome *p)
{
	register int i, j;
	register float fitness = 0;

	/* Fill evals */
	for (i=0; i<genelen; ++i) {
		evals[i].index = i;
		evals[i].key = p->prio[i];
	}

	/* Sort them (dumb sort for now) */
	for (i=0; i<genelen; ++i) {
		register float k = evals[i].key;

		for (j=i+1; j<genelen; ++j) {
			register float k2 = evals[j].key;

			if (k2 > k) {
				/* Swap */
				eval_t t = evals[i];
				evals[i] = evals[j];
				evals[j] = t;
				k = k2;
			}
		}
	}

	/* Fill-in the colors */
	for (i=0; i<genelen; ++i) {
		j = evals[i].index;
		pixels[(xdim * (point[j] / xdim)) + (point[j] % xdim)] =
		color[ evals[i].index ] = 1 + ((i * colors) / genelen);
	}

	/* Compute the fitness */
	for (i=0; i<genelen; ++i) {
		register int c = color[i];
		for (j=0; j<genelen; ++j) {
			if ((i != j) && (c == color[j])) {
				fitness += weight[i][j];
			}
		}
	}
	return(fitness);
}
	
void
write_P5(register geneome *p)
{
	/* Eval and output to stdout */

	printf("P5\n"
	       "# fitness is %f\n"
	       "%d %d\n"
	       "%d\n",
	       eval(p),
	       xdim,
	       ydim,
	       colors);
	fwrite(&(pixels[0]), 1, (xdim * ydim), stdout);
	fflush(stdout);
}

static inline float
randfloat(void)
{
	return((rand() & 0xffffff) / ((float) 0xffffff));
}

static void
randgene(register geneome *g)
{
	register int i;

	for (i=0; i<MAXLEN; ++i) {
		g->prio[i] = randfloat();
	}
}

int
biasrand(register int limit)
{
	return(rand() % (1 + (rand() % limit)));
}

static void
mutate(register geneome *kid,
register geneome *mom)
{
	/* Simple mutation */
	register int i;
	register int j = rand() % genelen;

	*kid = *mom;
	kid->prio[j] = randfloat();
}

static void
crossover(register geneome *kid,
register geneome *mom,
register geneome *dad)
{
	/* Simple crossover */
	register int i;

	for (i=0; i<genelen; ++i) {
		kid->prio[i] = ((rand() & 0x1000) ?
				mom->prio[i] :
				dad->prio[i]);
	}
}

main(int argc,
char **argv)
{
	register int i, j, k;
	register float v;
	register int pop, cross, mut, gens, gen = 0, dups = 0;
	register int whobest;
	register float vbest;

	if (argc == 6) {
		colors = atoi(argv[1]);
		pop = atoi(argv[2]);
		cross = atoi(argv[3]);
		mut = atoi(argv[4]);
		gens = atoi(argv[5]);
	} else {
		fprintf(stderr,
			"Usage: %s [colors] [pop] [cross] [mut] [gens]\n",
			argv[0]);
		exit(1);
	}

	srand(time(0));

	/* Read image */
	read_P5();

	/* The genetic algorithm loop... */
	do {
		register int a;

		/* Make a new guy... */
		if (gen < pop) {
			/* Still creating initial population */
			randgene( &(p[a = gen].g) );
		} else {
			/* Pick a pair of random pop members */
			register int b;

			a = (rand() % pop);
			do {
				b = (rand() % pop);
			} while (a == b);

			/* Make a be the worst one */
			if (p[a].v < p[b].v) {
				a ^= b; b ^= a; a ^= b;
			}

			/* Crossover or mutation? */
			if ((rand() % (cross + mut)) < cross) {
				/* Crossover needs a third */
				register int c = b;

				do {
					b = (rand() % pop);
				} while ((a == b) || (c == b));

				/* Make a be the worst one */
				if (p[a].v < p[b].v) {
					a ^= b; b ^= a; a ^= b;
				}

				/* Crossover */
				crossover(&(p[a].g),
					  &(p[b].g),
					  &(p[c].g));
			} else {
				/* Mutation */
				mutate(&(p[a].g), &(p[b].g));
			}
		}

		/* Evaluate the metric & update best */
		p[a].v = eval( &(p[a].g) );
		if ((gen == 0) || (p[a].v < vbest)) {
			whobest = a;
			vbest = p[a].v;

			/* Print the stats */
			fprintf(stderr,
				"Generation %5d: p[%d].v = %f\n",
				gen,
				whobest,
				vbest);
		}
	} while (++gen < gens);

	write_P5( &(p[whobest].g) );
}