/* ---------------------------------------------------------------------
    Celluvolve.
    Allan Crossman, Edinburgh iGEM 2011.

    This program is free software; you can redistribute it and/or
    modify it under the terms of the GNU General Public License
    as published by the Free Software Foundation; either version 2
    of the License, or (at your option) any later version.
--------------------------------------------------------------------- */
/* ---------------------------------------------------------------------
This program is a graphic demonstration of how cells with cellulases
attached to their cell surface will evolve high expression of said
cellulases faster than cells which simply release cellulases into the
media (in both cases assuming that producing more cellulase is useful).

The rules are as follows:

- An array of bacteria exists.
- Each bacterium has a "genotype": its ability to produce cellulase.
- The "genotype" is a number between 0 and 255.
- Each iteration, some random bacteria are chosen.
- These bacteria have a chance to reproduce, based on their fitness.
- Fitness is defined either as:
      * (LEFT) The ability of the cell to produce cellulase.
      * (RIGHT) The average ability of nearby cells (including itself).
- Reproduction means:
      * Overwriting a nearby cell with one's own genotype.
      * However, reproduction is not exact. There is "mutation".
      * The offspring's genotype is a number close to the parent's.

The left half of the simulation models the case where the cellulase
produced by a cell can only benefit itself. The right half models the
case where the cellulase can float away and benefit neighbours instead
of the cell that produced it.
--------------------------------------------------------------------- */

#include <assert.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <malloc.h>
#include <math.h>
#include <time.h>

/* Some constants that can reasonably be changed...

   RANDOMDISPERSAL - define this to just send the daughter to a random spot anywhere
   DISPERSALRADIUS - how far away a daughter can be placed
   NEIGHBOURSTOCONSIDER - probably doesn't matter too much
   CONSIDERATIONRADIUS - essentially, how far cellulases can travel. Lowering this
         helps the bacteria that produce free-floating cellulases.
*/

#undef RANDOMDISPERSAL
#define DISPERSALRADIUS 2	
#define NEIGHBOURSTOCONSIDER 5
#define CONSIDERATIONRADIUS 1

#define WIDTH 300
#define HEIGHT 300
#define SELECTIONS 1000

// How long to run the simulation?
#define FINALGENERATION 10000

// How often to report the average genotype?
#define REPORTTICK 50

// Should we draw BMP files when we report?
#define DRAWBMP

// Leave this alone...
#define DOUBLEWIDTH ((WIDTH * 2))

// The main data arrays...
int ability_attached[WIDTH][HEIGHT];
int ability_floating[WIDTH][HEIGHT];


float rnd (void)    // Random number between 0 and 1.
{
	return (float)rand() / RAND_MAX;
}

int intrnd (int max)    // Return integer between 0 and max inclusive.
{
	int result;

	result = rnd() * (max + 1);
	result %= max + 1;
	return result;
}

int getfitness (int x, int y)
{
	int sum = 0;
	int n = 0;
	int tarx; int tary;
	while (n < NEIGHBOURSTOCONSIDER)
	{
		tarx = x + intrnd(CONSIDERATIONRADIUS * 2) - CONSIDERATIONRADIUS;
		tary = y + intrnd(CONSIDERATIONRADIUS * 2) - CONSIDERATIONRADIUS;
		if (tarx >= 0 && tarx <= WIDTH - 1 && tary >= 0 && tary <= HEIGHT - 1)
		{
			n++;
			sum += ability_floating[tarx][tary];
		}
	}
	return sum / NEIGHBOURSTOCONSIDER;
}

float getaverage_attached (void)
{
	int sum = 0;
	int x = 0;
	int y = 0;

	for (x = 0; x < WIDTH; x++)
	{
		for (y = 0; y < HEIGHT; y++)
		{
			sum += ability_attached[x][y];
		}
	}
	return (float) sum / (WIDTH * HEIGHT);
}

float getaverage_floating (void)
{
	int sum = 0;
	int x = 0;
	int y = 0;

	for (x = 0; x < WIDTH; x++)
	{
		for (y = 0; y < HEIGHT; y++)
		{
			sum += ability_floating[x][y];
		}
	}
	return (float) sum / (WIDTH * HEIGHT);
}

void drawbmp (char * filename)
{ 
	unsigned int headers[13];
	FILE * outfile;
	int extrabytes;
	int paddedsize;
	int x; int y; int n;
	int red, green, blue;
 
	extrabytes = 4 - ((DOUBLEWIDTH * 3) % 4);		// How many bytes of padding to add to each
									// horizontal line - the size of which must
									// be a multiple of 4 bytes.
	if (extrabytes == 4) extrabytes = 0;
 
	paddedsize = ((DOUBLEWIDTH * 3) + extrabytes) * HEIGHT;
 
	// Headers...
	// Note that the "BM" identifier in bytes 0 and 1 is NOT included in these "headers".
 
	headers[0]  = paddedsize + 54;	// bfSize (whole file size)
	headers[1]  = 0;				// bfReserved (both)
	headers[2]  = 54;				// bfOffbits
	headers[3]  = 40;				// biSize
	headers[4]  = DOUBLEWIDTH;		// biWidth
	headers[5]  = HEIGHT;			// biHeight
 
	// Would have biPlanes and biBitCount in position 6, but they're shorts.
	// It's easier to write them out separately (see below) than pretend
	// they're a single int, especially with endian issues...
 
	headers[7]  = 0;				// biCompression
	headers[8]  = paddedsize;		// biSizeImage
	headers[9]  = 0;				// biXPelsPerMeter
	headers[10] = 0;				// biYPelsPerMeter
	headers[11] = 0;				// biClrUsed
	headers[12] = 0;				// biClrImportant
 
	outfile = fopen(filename, "wb");
 
	//
	// Headers begin...
	// When printing ints and shorts, we write out 1 character at a time to avoid endian issues.
	//
 
	fprintf(outfile, "BM");
 
	for (n = 0; n <= 5; n++)
	{
		fprintf(outfile, "%c", headers[n] & 0x000000FF);
		fprintf(outfile, "%c", (headers[n] & 0x0000FF00) >> 8);
		fprintf(outfile, "%c", (headers[n] & 0x00FF0000) >> 16);
		fprintf(outfile, "%c", (headers[n] & (unsigned int) 0xFF000000) >> 24);
	}
 
	// These next 4 characters are for the biPlanes and biBitCount fields.
 
	fprintf(outfile, "%c", 1);
	fprintf(outfile, "%c", 0);
	fprintf(outfile, "%c", 24);
	fprintf(outfile, "%c", 0);
 
	for (n = 7; n <= 12; n++)
	{
		fprintf(outfile, "%c", headers[n] & 0x000000FF);
		fprintf(outfile, "%c", (headers[n] & 0x0000FF00) >> 8);
		fprintf(outfile, "%c", (headers[n] & 0x00FF0000) >> 16);
		fprintf(outfile, "%c", (headers[n] & (unsigned int) 0xFF000000) >> 24);
	}
 
	//
	// Headers done, now write the data...
	//
 
	for (y = HEIGHT - 1; y >= 0; y--)		// BMP image format is written from bottom to top...
	{
		for (x = 0; x <= DOUBLEWIDTH - 1; x++)
		{
 			if (x < WIDTH)
			{
				red = ability_attached[x][y];
				green = ability_attached[x][y];
				blue = ability_attached[x][y] / 2;
			} else {
				red = ability_floating[x - WIDTH][y];
				green = ability_floating[x - WIDTH][y];
				blue = ability_floating[x - WIDTH][y] / 2;
			}
 
			// Also, it's written in (b,g,r) format...
 
			fprintf(outfile, "%c", blue);
			fprintf(outfile, "%c", green);
			fprintf(outfile, "%c", red);
		}
		if (extrabytes)		// See above - BMP lines must be of lengths divisible by 4.
		{
			for (n = 1; n <= extrabytes; n++)
			{
				fprintf(outfile, "%c", 0);
			}
		}
	}
 
	fclose(outfile);
	return;
}

int main (void)
{
	int x; int y;
	int tarx; int tary;
	int n;
	int fitness;
	int tick = 0;
	char filename[200];

	printf("\nCelluvolve, an iGEM model from Edinburgh 2011...\n\n");

	for (x = 0; x < WIDTH; x++)
	{
		for (y = 0; y < HEIGHT; y++)
		{
			ability_attached[x][y] = 50;
			ability_floating[x][y] = 50;
		}
	}

	srand((int)time(NULL));

	printf("Iteration\tAverage genotype (left)\tAverage genotype (right)\n");
	while (tick < FINALGENERATION)
	{
		tick++;

		// Fitness is your own ability...
		for (n = 0; n < SELECTIONS; n++)
		{
			x = intrnd(WIDTH - 1);
			y = intrnd(HEIGHT - 1);
			fitness = ability_attached[x][y];
			if (intrnd(255) < fitness)
			{
				#ifdef RANDOMDISPERSAL
					tarx = intrnd(WIDTH - 1);
					tary = intrnd(HEIGHT - 1);
				#else
					tarx = x + intrnd(DISPERSALRADIUS * 2) - DISPERSALRADIUS;
					tary = y + intrnd(DISPERSALRADIUS * 2) - DISPERSALRADIUS;
				#endif
				if (tarx >= 0 && tarx <= WIDTH - 1 && tary >= 0 && tary <= HEIGHT - 1)
				{
					ability_attached[tarx][tary] = ability_attached[x][y] + intrnd(24) - 12;
					if (ability_attached[tarx][tary] < 0) ability_attached[tarx][tary] = 0;
					if (ability_attached[tarx][tary] > 255) ability_attached[tarx][tary] = 255;
				}
			}
		}
		// Fitness is the average ability of you and your neighbours...
		for (n = 0; n < SELECTIONS; n++)
		{
			x = intrnd(WIDTH - 1);
			y = intrnd(HEIGHT - 1);
			fitness = getfitness(x, y);
			if (intrnd(255) < fitness)
			{
				#ifdef RANDOMDISPERSAL
					tarx = intrnd(WIDTH - 1);
					tary = intrnd(HEIGHT - 1);
				#else
					tarx = x + intrnd(DISPERSALRADIUS * 2) - DISPERSALRADIUS;
					tary = y + intrnd(DISPERSALRADIUS * 2) - DISPERSALRADIUS;
				#endif
				if (tarx >= 0 && tarx <= WIDTH - 1 && tary >= 0 && tary <= HEIGHT - 1)
				{
					ability_floating[tarx][tary] = ability_floating[x][y] + intrnd(24) - 12;
					if (ability_floating[tarx][tary] < 0) ability_floating[tarx][tary] = 0;
					if (ability_floating[tarx][tary] > 255) ability_floating[tarx][tary] = 255;
				}
			}
		}
		if (tick % REPORTTICK == 0)
		{
			printf("%6d\t%5f\t%5f\n", tick, getaverage_attached(), getaverage_floating());
			#ifdef DRAWBMP
				sprintf(filename, "%05d.bmp", tick);
				drawbmp(filename);
			#endif
		}
	}

	return 0;
}