/**
 * @file ga.c
 * @brief Genetic Algorithm Core
 * @author Telford Tendys
 *
 * @date 2012-06-10
 *
 * Simulation of people trading with each other, based on the presumption of a
 * "Prisoner's Dilemma" in a continuous space (makes a sort of saddle shape).
 * Thus either player can cooperate more (maximum 1.0) or less (minimum 0.0)
 * and thus we get payoffs for the encounter. In addition, players have a bitmap
 * that allows them to compare themselves to other players and a tolerance that
 * acts as a bias (players with high tolerance are likely to cooperate with other
 * players even those quite different to themselves, players with low tolerance
 * will only cooperate with players very similar to themselves).
 *
 * Copyright (C) 2010, 2011 Telford Tendys <telford@lnx-bsp.net>
 *
 * This code is distributed under the Creative Commons "Attribution-ShareAlike 3.0"
 * SEE: http://creativecommons.org/licenses/by-sa/3.0/
 *
 * At your option you may convert this to the Free Software Foundation's GNU
 * General Public License version 3 or any later version.
 * SEE: http://www.gnu.org/licenses/gpl.html
 */

/** @mainpage Quickstart to GA Internals
 *
 * @section intro Introduction
 *
 * <a href=http://evolve.lnx-bsp.net/>This article outlines the fundamental process.</a>
 *
 * As a simple outline, all of the program is in one single file. There are various
 * main() functions for test modes (don't worry about them right now) but the actual main()
 * requires -DRUN_MODE compile flag. This function does the following:
 * <ul>
 * <li>Set up a new cellular grid.
 * <li>Iterate a fixed number of cycles.
 * <li>Pick pairs of nearby cells.
 * <li>Evaluate an interaction between pairs.
 * <li>Save frames as PNM images every so often.
 * </ul>
 *
 * That's about the heart of it.
 *
 * @section index Key Algorithm Functions and Structures
 *
 * <ul>
 * <li>struct cell_S contains the data held in each cell.
 * <li>function point_select() picks a random pair of nearby cells.
 * <li>function hamming32() calculates Hamming distance.
 * <li>function tolfunc() calculates the parametric tolerance (based on Hamming distance).
 * <li>function outcome() calculates the outcome of the Prisoner's Dilemma.
 * <li>function transaction() is the full interaction between the cells.
 * </ul>
 *
 * @section index Image and Colour Related
 *
 * <ul>
 * <li>struct col_item_S contains a single item mapping from Hamming space to RGB colour space.
 * <li>array col_item_T col_items[] is a list of standard colours and bitmaps.
 * <li>function get_colour() converts from Hamming space to RGB colour space (using col_items).
 * <li>function output_frame() dumps the whole cell grid to an image file (with colour translation, etc).
 * </ul>
 *
 */

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

/** Total width of the cellular grid (also output image width in pixels) */
#define WIDTH 640
/** Total height of the cellular grid (also output image height in pixels) */
#define HEIGHT 480

/** The penalty (i.e. a negative payoff) in return for offering best attitude (1.0) in response to worst attitude (0.0) */
const double PAYOFF_SUCKER       = -20;

/** The penalty (i.e. a negative payoff) in return for both sides offering the worst attitude (0.0) */
const double PAYOFF_WAR          = -15;

/** The penalty (i.e. a negative payoff) in return for both sides offering the worst attitude (1.0) */
const double PAYOFF_PEACE        =   5;

/** The penalty (i.e. a negative payoff) in return for offering worst attitude (0.0) in response to best attitude (1.0) */
const double PAYOFF_TEMPTATION   =  10;

/** Initial value for cell_S::min_split when system bootstraps */
const double INITIAL_MIN_SPLIT   =  100;

/** Regular maintenance cost for all cells, regardless of attitude */
const double MAINTENANCE         =  -3;

/** Initial value for cell_S::res when system bootstraps */
const double INITIAL_RES         =  2000;

/** Initial value for both cell_S::tol_low and cell_S::tol_high when system bootstraps */
const double INITIAL_TOL         =  0.75;

/** Initial value for cell_S::mobility when system bootstraps */
const double INITIAL_MOB         =  0.00;

/** Termination condition for the simulation, can be quite a large number */
const long long TOTAL_CYCLES     =  40000000000LL;

/** Artificial upper limit on the amount of resources any cell can  */
const double MAX_RES             =  5000;

/** Factor for random mistakes and confusion amongst cells.
 * A bigger value here means cells are more likely to misjudge each other when evaluating trust
 * (e.g. cultural patterns might not be properly recognised).
 */
const double CONFUSION_FACTOR_1  =  10;

/** Factor for random mistakes and confusion amongst cells.
 * A bigger value here means cells are more likely to misjudge their own reactions
 * (e.g. over react or under react to an encounter).
 */
const double CONFUSION_FACTOR_2  =  10;

/**
 * Artificial upper limit on high-mobility cells, possibly because of some physical
 * limitation or whatever. If set to 0.0 then all mobility is impossible, but
 * the highest value is 1.0 and cells are allowed to select for LOWER mobility
 * than this value (although experimentally, it seems to always result in a spread
 * from zero up to the limit value).
 */
const double MOBILITY_BLOCKER = 0.3;

/** A frame counter for output images, not all frames are used in the final video.
    The frame counter goes into the filename */
int frame = 0;

/** A counter that controls the length of the simulation */
long long cycle_count = -1;

/* -------------------------------------------------------------------------------- */
/*                                 random generator                                 */
/* -------------------------------------------------------------------------------- */

typedef int TS32;                   /* Presumes 32 bit integers -- the gcc standard */

/** Maximum size of random state data */
#define TEL_RAND_SIZE 6

/** Random generator static state data */
typedef TS32 TRND32[ TEL_RAND_SIZE ];

/** Used by random generator internals */
TS32 *base_hash_matrix = 0;

static int prep_hash_base( void )
{
	TS32 i, j, k;

	base_hash_matrix = malloc( 4 * 256 );
	if( !base_hash_matrix )
	{
		fprintf( stderr, "Out of memory!\n" );
		exit( -1 );
	}
	bzero( base_hash_matrix, 4 * 256 );
	j = 0; k = 0;
	for( i = 0; i < 100000; i++ )
	{
		j += i;                             // Triangular numbers
		k += j + i;                         // Tetrahedral numbers
		base_hash_matrix[ i & 0xFF ] += k;  // Hyper-tetrahedral numbers
	}
	return( 0 );
}

/** Extract the low byte of an integer (0 to 255) as unsigned */
#define LOW_BYTE(x) ((unsigned char)(x))

/**
 * Generate a random value, based on the given array of state
 * variables (which are also updated by the random generation).
 * Faster than the glibc random generator, possibly not quite as
 * random but pretty good all told.
 *
 * Output value is in rand_state[ 0 ]
 *
 * @return always 0
 */
int tel_random3( TRND32 rand_state )
{
	int e;
	register TS32 state;

	if( !base_hash_matrix ) { e = prep_hash_base(); if( e < 0 ) return( e ); }

// ----------------------------------------
	state = rand_state[ 1 ];
	state >>= 1;
	state ^= base_hash_matrix[ LOW_BYTE( state )];
	rand_state[ 1 ] = state;
	state += rand_state[ 2 ];
	state >>= 1;
	state ^= base_hash_matrix[ LOW_BYTE( state )];
	rand_state[ 2 ] = state;
	state += rand_state[ 3 ];
	state >>= 1;
	state ^= base_hash_matrix[ LOW_BYTE( state )];
	rand_state[ 3 ] = state;
	state += rand_state[ 4 ];
	state >>= 1;
	state ^= base_hash_matrix[ LOW_BYTE( state )];
	rand_state[ 4 ] = state;
// ----------------------------------------
	state >>= 8;
	state ^= rand_state[ 3 ];
	state >>= 8;
	state ^= rand_state[ 2 ];
	state >>= 8;
	state ^= rand_state[ 1 ];
	rand_state[ 0 ] = state;

	return( 0 );
}

/** Static random state variable, starts all zero when program begins */
TRND32 rand_state;
/**
 * Use the tel_random3() function to supply a stream of
 * pseudo-random numbers, using a static state variable.
 *
 * @return pseudo-random 32 bit number (signed)
 */
TS32 rand32( void )
{
	tel_random3( rand_state );
	return( rand_state[ 0 ]);
}

/**
 * Take a probability value (forced to limit between 0 and 1)
 * and calculate whether an event happens (return 1 if it happens 0 otherwize)
 */
int random_chance( double p )
{
	if( p <= 0 ) { return( 0 ); }
	if( p >= 1 ) { return( 1 ); }
	int x = rand32();
	if( x < 0 ) { x = -x; }
	p *= 0x7FFFFFFF;
	if( p > x ) { return( 1 ); }
	return( 0 );
}

/* -------------------------------------------------------------------------------- */
/*                                 hamming distance                                 */
/* -------------------------------------------------------------------------------- */

/**
 * A lookup table for calculating the Hamming distance between two bytes.
 * Counts the number of 1 bits in a single byte.
 */
unsigned char byte_dis[] =
{
    0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4,
    1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5,
    1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5,
    2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
    1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5,
    2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
    2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
    3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
    1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5,
    2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
    2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
    3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
    2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
    3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
    3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
    4, 5, 5, 6, 5, 6, 6, 7, 5, 6, 6, 7, 6, 7, 7, 8,
};

/**
 * This is merely a typedef to make it easier to type "struct cell_S"
 * but it means nothing more than that.
 */
typedef struct cell_S cell_T;

/**
 * The cell structure is the heart of the evolutionary simulation
 * because these are the parameters which can evolve, and this data
 * is stored in the cellular grid.
 */
struct cell_S
{
	/** Total available resources in this particular cell */
	double res;

	/** Tolerance function parameter (low value) -- see tolfunc() */
	double tol_low;

	/** Tolerance function parameter (mid value) -- see tolfunc() */
	double tol_mid;

	/** Tolerance function parameter (high value) -- see tolfunc() */
	double tol_high;

	/** Estimate of local region diversity (based on interractions with this cell) */
	double diversity;

	/** Threshold for splitting the cell (asexual reproduction) which requires resources and
		also available empty space to reproduce into. */
	double min_split;

	/** Probability of cell swapping places with neighbour during transaction. */
	double mobility;

	/** Hamming-space bit pattern representing this cell's "cultural itentity" */
	unsigned int bits;
};

/**
 * Evaluate hamming distance between two 32bit numbers
 *
 * @return a number between 0 and 32
 */
unsigned int hamming32( unsigned int x1, unsigned int x2 )
{
	unsigned int tot;

	x1 ^= x2;
	tot = byte_dis[ x1 & 0xFF ];
	x1 >>= 8;
	tot += byte_dis[ x1 & 0xFF ];
	x1 >>= 8;
	tot += byte_dis[ x1 & 0xFF ];
	x1 >>= 8;
	tot += byte_dis[ x1 & 0xFF ];
	return( tot );
}

/* -------------------------------------------------------------------------------- */
/*                               tolerance calculation                              */
/* -------------------------------------------------------------------------------- */

/**
 * Parametric parabola made with shape functions and trapped between 0, 1
 *
 * Allows significant scope for the tolerance function, including high-side, low-side and
 * how much curve there is in the middle. This allows the agents to tune their attitude to
 * other agents, including some degree of non-linearity.
 */
double tolfunc( double x, double q_l, double q_m, double q_h )
{
	double y = 0;

	if( x < 0 ) x = 0;
	if( x > 1 ) x = 1;
	y +=     ( 1 - x ) * q_l; /* Param #1 adjusts low side  */
	y += x * ( 1 - x ) * q_m; /* Param #2 adjusts middle    */
	y += x             * q_h; /* Param #3 adjusts high side */
	if( y > 1 ) y = 1;
	if( y < 0 ) y = 0;
        return( y );
}

/**
 * A random pruturbation (or mutation if you prefer) that allows the evolutionary
 * algorithm to continue
 */
double adjustment( void )
{
	double x = rand32() - rand32();
	x *= 1e-11;
//	printf( "%f\n", x );
	return( x );
}

/**
 * return [ 0.0 to 1.0 ] where 0.0 is the worst and 1.0 is the best. tol values is taken
 * from A point of view, and the above tolfunc().
 */
double att( cell_T *A, cell_T *B )
{
	double ham = hamming32( A->bits, B->bits );
	double ret;

	ham /= 32;
	ham += adjustment() * CONFUSION_FACTOR_1;
	ret = tolfunc( ham, A->tol_low, A->tol_mid, A->tol_high );
#ifdef RUN_MODE
	/*
	 * We output periodic samples of this calculation for statistical analysis later
	 */
	if(( cycle_count & 0x0FFFF ) == 0 )
	{
		printf( "%lld\t%f\t%f\t%f\t%f\t%f\t%f\t%f\n", cycle_count, ham, A->tol_low, A->tol_mid, A->tol_high, ret, A->mobility, A->res );
	}
#endif
	return( ret + adjustment() * CONFUSION_FACTOR_2 );
}

#ifdef TEST_MODE_1

int main()
{
	/*
	 * Testing the att system based on artifical case
	 */
	cell_T A = {0};
	cell_T B = {0};

	A.tol_low = 0;
	A.tol_high = 1;

	for( A.tol_mid = -2; A.tol_mid <= 2; A.tol_mid += 0.1 )
	{
		for( A.bits = 0xFFFFFFFF; ; A.bits >>= 1 )
		{
			printf( "%f\t%f\t%f\n", A.tol_mid, 1.0 - ( hamming32( A.bits, B.bits ) / 32.0 ), att( &A, &B ));
			if( A.bits == 0 ) break;
		}
	}
	return( 0 );
}

#endif // TEST_MODE_1

/* -------------------------------------------------------------------------------- */
/*                            interaction outcome calculation                       */
/* -------------------------------------------------------------------------------- */

/**
 * Evaluate the outcome of the Prisoner's Dilemma game by taking two attitudes
 * and linearly interpolating between the four extreme possible outcomes:
 * PAYOFF_PEACE, PAYOFF_TEMPTATION, PAYOFF_SUCKER, PAYOFF_WAR
 *
 * All inputs are bounded between 0 and 1 (inclusive) so callers to this function
 * need not be overly cautious about checking their values are within range.
 */
double outcome( double att_A, double att_B )
{
	double z;

	if( att_A < 0 ) { att_A = 0; }
	if( att_A > 1 ) { att_A = 1; }
	if( att_B < 0 ) { att_B = 0; }
	if( att_B > 1 ) { att_B = 1; }

	z = ( 0
		+     att_A     *     att_B     * PAYOFF_PEACE      /* Both parties are helpful  */
		+ ( 1 - att_A ) *     att_B     * PAYOFF_TEMPTATION /* Hostility vs helpful      */
		+     att_A     * ( 1 - att_B ) * PAYOFF_SUCKER     /* Helpful vs hostile        */
		+ ( 1 - att_A ) * ( 1 - att_B ) * PAYOFF_WAR        /* Both parties are hostile  */
	);

	return( z );
}

#ifdef TEST_MODE_1A

/*
 *	Scan the whole outcome() function.
 *
 *	This makes it good for a plot, etc.
 */

int main()
{
	double x, y;

	printf( "\"x\"\t\"y\"\t\"v\"\n" );
	for( x = 0; x <= 1.005; x += 0.01 )
	{
		if( x > 1 ) x = 1;
		for( y = 0; y <= 1.005; y += 0.01 )
		{
			if( y > 1 ) y = 1;
			printf( "%f\t%f\t%f\n", x, y, outcome( x, y ));
		}
	}
	return( 0 );
}

#endif // TEST_MODE_1A


#ifdef TEST_MODE_1B

/*
 * Output hamming distance from random pairs
 */
int main()
{
	unsigned int A, B, i;
	for( i = 0; i < 1000000; i++ )
	{
		A = rand32();
		B = rand32();
		printf( "%d\n", hamming32( A, B ));
	}
	return( 0 );
}

#endif // TEST_MODE_1B

/**
 * Split one cell into two.
 *
 * All inheritable traits are copied exactly.
 * Resources are split in half and each side gets the same.
 * Cell division is an instantaneous operation so there is no stage where the cell
 * is "out of operation" due to the reproductive pricess.
 */
void do_split( cell_T *A, cell_T *B )
{
	A->res /= 2;
	B->res = A->res;
	B->tol_low  = A->tol_low;
	B->tol_mid  = A->tol_mid;
	B->tol_high = A->tol_high;
	B->bits = A->bits;
	B->min_split = A->min_split;
	B->mobility = A->mobility;
}

/**
 * Evaluate one transaction between a pair of cells
 */
void transaction( cell_T *A, cell_T *B )
{
	double att_AB = att( A, B );
	double att_BA = att( B, A );
	double return_A = 0;
	double return_B = 0;

	A->res += MAINTENANCE;

	/* Game outcome -- e.g. prisoner's dilemma */
	if( A->res > 0 && B->res > 0 )
	{
		return_A = outcome( att_AB, att_BA );
		return_B = outcome( att_BA, att_AB );
		A->res += return_A;
		B->res += return_B;
	}

	if( A->res > MAX_RES ) { A->res = MAX_RES; }
	if( B->res > MAX_RES ) { B->res = MAX_RES; }

	/* Not strictly necessary but keep it neat */
	if( A->res < 0 ) { A->res = 0; }
	if( B->res < 0 ) { B->res = 0; }

	/* Takeover, one cell splits to take a vacant space */
	if( A->res <= 0 && B->res > B->min_split ) do_split( B, A );
	if( B->res <= 0 && A->res > A->min_split ) do_split( A, B );

	/* Random mutation (A only, non-vacant only) */
	if( A->res <= 0 ) return;

#if 1
	{

		int mask = 0;
		switch( 0x0FFF & rand32())
		{
			case  0: mask = 0x00000001; break;
			case  1: mask = 0x00000002; break;
			case  2: mask = 0x00000004; break;
			case  3: mask = 0x00000008; break;

			case  4: mask = 0x00000010; break;
			case  5: mask = 0x00000020; break;
			case  6: mask = 0x00000040; break;
			case  7: mask = 0x00000080; break;

			case  8: mask = 0x00000100; break;
			case  9: mask = 0x00000200; break;
			case 10: mask = 0x00000400; break;
			case 11: mask = 0x00000800; break;

			case 12: mask = 0x00001000; break;
			case 13: mask = 0x00002000; break;
			case 14: mask = 0x00004000; break;
			case 15: mask = 0x00008000; break;

			case 16: mask = 0x00010000; break;
			case 17: mask = 0x00020000; break;
			case 18: mask = 0x00040000; break;
			case 19: mask = 0x00080000; break;

			case 20: mask = 0x00100000; break;
			case 21: mask = 0x00200000; break;
			case 22: mask = 0x00400000; break;
			case 23: mask = 0x00800000; break;

			case 24: mask = 0x01000000; break;
			case 25: mask = 0x02000000; break;
			case 26: mask = 0x04000000; break;
			case 27: mask = 0x08000000; break;

			case 28: mask = 0x10000000; break;
			case 29: mask = 0x20000000; break;
			case 30: mask = 0x40000000; break;
			case 31: mask = 0x80000000; break;
		}

		if( mask )
		{
			/* Some of the A->bits might change to either align or disalign with B->bits */
			if( return_A > 0.1 )
			{
				int bits = A->bits ^ B->bits;
				bits &= mask;
				A->bits ^= bits; // Align bit in A to match B
			}

			if( return_A < 0.1 )
			{
				int bits = ~( A->bits ^ B->bits );
				bits &= mask;
				A->bits ^= bits; // Disalign bit in A to be different to B
			}
		}
	}
#endif

	// Allow A the choice to swap places with B and thus put mobility into the system
	if( random_chance( A->mobility ))
	{
		// Swap the pair
		cell_T tmp;
		bcopy(    A, &tmp, sizeof( tmp ));
		bcopy(    B,    A, sizeof( tmp ));
		bcopy( &tmp,    B, sizeof( tmp ));
	}

	A->tol_low += adjustment();
	A->tol_mid += adjustment();
	A->tol_high += adjustment();
	A->min_split += adjustment();
	A->mobility += adjustment();
	if( A->min_split < 0 ) { A->min_split = 0; }
	// Since mobility is a probability value, less than zero doesn't even make sense.
	if( A->mobility < 0 ) { A->mobility = 0; }
	// We force an arbitrary upper limit, in order to discover what influence
	// mobility has on the overall system dynamics.
	if( A->mobility > MOBILITY_BLOCKER ) { A->mobility = MOBILITY_BLOCKER; }
}

/**
 * Load initial values into a cell at bootstrap time
 */
void setup( cell_T *A )
{
	A->res      = INITIAL_RES;
	A->tol_low  = INITIAL_TOL;
	A->tol_mid  = 0.00;         /* Start with flat curve */
	A->tol_high = INITIAL_TOL;
	A->bits = rand32();
	A->min_split = INITIAL_MIN_SPLIT;
	A->mobility = INITIAL_MOB;
}

/**
 * This is the cellular grid
 */
cell_T grid[ WIDTH * HEIGHT ];

/**
 * Pick a random pair of cells, that must be neighbours or near-neighbours.
 * More likely to pick close neighbours than distant neighbours.
 * FIXME: This could probably be better.
 */
void point_select( cell_T **Ap, cell_T **Bp )
{
	for(;;)
	{
		int j = ( rand32() & 0x7FFFF );
		int k = j;

		if( j >= ( WIDTH * HEIGHT )) { continue; }
		switch( rand32() % 12 )
		{
			case  0: k++;                break;
			case  1: k--;                break;
			case  2: k += WIDTH;         break;
			case  3: k -= WIDTH;         break;
			case  4: k++; k += WIDTH;    break;
			case  5: k--; k += WIDTH;    break;
			case  6: k++; k -= WIDTH;    break;
			case  7: k--; k -= WIDTH;    break;

			case  8: k += 2;             break;
			case  9: k -= 2;             break;
			case 10: k += 2*WIDTH;       break;
			case 11: k -= 2*WIDTH;       break;
		}
		if( k < 0 || k >= ( WIDTH * HEIGHT )) { continue; }
		*Ap = grid + j;
		*Bp = grid + k;
		return;
	}
}

/* -------------------------------------------------------------------------------- */

#ifdef TEST_MODE_2

int main()
{
	/*
	 * Testing the att system based on artifical case
	 */
	int i;
	cell_T A = {0};
	cell_T B = {0};

	rand_state[ 1 ] = time(0);

	setup( &A );
	setup( &B );

	for( i = 0; i < 1000000; i++ )
	{
		transaction( &A, &B );
		transaction( &B, &A );
		printf( "%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\n",
			A.res, A.tol_low, A.tol_mid, A.tol_high,
			B.res, B.tol_low, B.tol_mid, B.tol_high,
			hamming32( A.bits, B.bits ) / 32.0 );
	}
	return( 0 );
}

#endif // TEST_MODE_2


#ifdef TEST_MODE_3

int main()
{
	printf( "%f\t%f\t%f\t%f\n",
		outcome( 0, 0 ),
		outcome( 0, 1 ),
		outcome( 1, 0 ),
		outcome( 1, 1 ));
	return( 0 );
}

#endif


#ifdef TEST_MODE_4

int main()
{
	/*
	 * Testing the random selection of points
	 */
	int i;
	int frame = 0;

//	rand_state[ 1 ] = time(0);

	for( i = 0; i < 1000000000; i++ )
	{
		cell_T *A;
		cell_T *B;

		point_select( &A, &B );
		A->res++;
		B->res--;

		if(( i & 0x03FFFFF ) == 0 )
		{
			FILE *fp;
			char buf[ 256 ];
			int j, k;

			sprintf( buf, "pic_%04d.ppm", frame++ );
			fp = fopen( buf, "w" );
			fprintf( fp, "P6\n%d %d\n255\n", WIDTH, HEIGHT );

			A = grid;
			for( j = 0; j < HEIGHT; j++ )
			{
				for( k = 0; k < WIDTH; k++ )
				{
					double r = 0, g = 0, b = 0;

					r = A->res;                // Relative visits A - B
					r /= 2;
					r = floor( r + 127.5 );
					if( r < 0 ) r = 0;
					if( r > 255 ) r = 255;

					fputc(( int )r, fp );
					fputc(( int )g, fp );
					fputc(( int )b, fp );
					A++;
				}
			}
			fclose( fp );
			fprintf( stderr, "frame=%d\n", frame ); fflush( stdout );
		}
	}
	return( 0 );
}

#endif // TEST_MODE_4


/* -------------------------------------------------------------------------------- */
/*                            output colour calculation                             */
/* -------------------------------------------------------------------------------- */

/**
 * This is merely a typedef to make it easier to type "struct col_item_S"
 * but it means nothing more than that.
 */
typedef struct col_item_S col_item_T;

/**
 * One colour item, mapping from a bit pattern to some ( r, g, b ) tuplet.
 * Note that r, g and b are doubles with values from 0 to 1 (inclusive).
 */
struct col_item_S
{
	/** The Hamming space bit vector representing this particular colour */
	unsigned int bits;
	/** The (r, g, b) tuplet representing the colour in RGB space */
	float r, g, b;
};

/**
 * A list of pre-defined colours are here.
 *
 * These are largely arbitrary but they span a range of primary and secondary RBG
 * colour values (i.e. bright coulours) and they span various orthogonal bitmaps.
 */
col_item_T col_items[] = {
	{ 0x00000000, 1, 0, 0 },
	{ 0xFFFFFFFF, 0, 1, 1 },
	{ 0xFFFF0000, 0, 1, 0 },
	{ 0x0000FFFF, 1, 0, 1 },
	{ 0x00FFFF00, 0, 0, 1 },
	{ 0xFF0000FF, 1, 1, 0 },
	{ 0x3C3C3C3C, 1, 1, 1 },
};

/** The number of colours in col_items */
#define N_COL 7

/**
 * Convert bit pattern to RGB colour for display purposes.
 *
 * This does a fairly simple merge of multiple vectors using the best and second-best
 * fit against a fixed list of vectors (based on searching for the shortest Hamming distance).
 *
 * The output R, G, B values determines the colour for this particular cell,
 * although colours are largely arbitrary, they make it easy to visualise the boundary of
 * cooperative groups, and also tend to reveal the overall diversity within the group.
 * Of course, collapsing a 32 dimensional space down to a 3 dimensional space will
 * necessarily lose some information.
 */

int get_colour( unsigned int bits, double *r, double *g, double *b )
{
	unsigned int j;
	int d1, d2, j1, j2;
	double ratio;

	j1 = j2 = -1;
	d1 = d2 = 100;
	for( j = 0; j < N_COL; j++ )
	{
		int dis = hamming32( bits, col_items[ j ].bits );
		if( dis < d1 )
		{
			d2 = d1;
			j2 = j1;
			d1 = dis;
			j1 = j;
		}
		else if( dis < d2 )
		{
			d2 = dis;
			j2 = j;
		}
	}
	ratio = d1;
	ratio /= ( d1 + d2 );
	*r = col_items[ j1 ].r * ( 1 - ratio ) + col_items[ j2 ].r * ratio;
	*g = col_items[ j1 ].g * ( 1 - ratio ) + col_items[ j2 ].g * ratio;
	*b = col_items[ j1 ].b * ( 1 - ratio ) + col_items[ j2 ].b * ratio;

#ifdef TEST_MODE_5
	printf( "d1=%d j1=%d d2=%d j2=%d\n", d1, j1, d2, j2 );
	printf( "ratio = %f\n", ratio );
	printf( "rgb = %.03f %.03f %.03f\n", *r, *g, *b );
#endif

	return( 0 );
}

#ifdef TEST_MODE_5

int main( int argc, char *argv[])
{
	int i;
	for( i = 1; i < argc; i++ )
	{
		unsigned int bits = strtoll( argv[ i ], 0, 0 );
		double r, b, g;

		get_colour( bits, &r, &g, &b );
		printf( "rgb = %.03f %.03f %.03f\n", r, g, b );
	}
	return( 0 );
}

#endif // TEST_MODE_5

#ifdef TEST_MODE_6

/*
 * This is just to make plots of the tolfunc() as based on a selection
 * of parameter values.
 */
double q_l, q_m, q_h;

int scan_x()
{
	double x;

	for( x = 0; x <= 1.001; x += 0.01 )
	{
		double y = tolfunc( x, q_l, q_m, q_h );
		printf( "%f\t%f\t%f\t%f\t%f\n", q_l, q_m, q_h, x, y );
	}
	// Separator record, to break the plot lines
	printf( "NA\tNA\tNA\tNA\tNA\n" );
}

int main( int argc, char *argv[])
{
	q_l = 0.5;
	q_m = 0.0;
	q_h = 0.5;

	for( q_l = -0.2; q_l <= 1.201; q_l += 0.2 ) { scan_x(); }
	q_l = 0.5;
	for( q_m = -2.0; q_m <= 2.01; q_m += 0.5 ) { scan_x(); }
	q_m = 0.0;
	for( q_h = -0.2; q_h <= 1.201; q_h += 0.2 ) { scan_x(); }
	return( 0 );
}

#endif // TEST_MODE_6

/* -------------------------------------------------------------------------------- */


#ifdef TEST_MODE_7
/*
 * List of the typical confusion adjustments, for analysis
 */

int main( int argc, char *argv[])
{
	int i;
	for( i = 0; i < 10000; i++ )
	{
		printf( "%f\t%f\n", adjustment() * CONFUSION_FACTOR_1, adjustment() * CONFUSION_FACTOR_2 );
	}
}

#endif // TEST_MODE_7

/* -------------------------------------------------------------------------------- */

#ifdef RUN_MODE

/**
 * Output a frame in PPM format which we can convert to something more practical later
 * (e.g. MPEG video stream). Colour determines approximate allegiance to some group.
 * Overall output intensity determies resource levels (black for dead cells, dark for near-dead).
 *
 * The PPM format is uncompressed and inefficient, but relatively easy to generate.
 */
int output_frame( void )
{
	FILE *fp;
	char buf[ 256 ];
	int j, k;
	cell_T *A = grid;

	sprintf( buf, "pic_%05d.ppm", frame++ );
	fp = fopen( buf, "w" );
	fprintf( fp, "P6\n%d %d\n255\n", WIDTH, HEIGHT );
	for( j = 0; j < HEIGHT; j++ )
	{
		for( k = 0; k < WIDTH; k++ )
		{
			double r = 0;
			double g = 0;
			double b = 0;
			double res;

			res = A->res;
			if( res > 0 )
			{
				// Try to do something of a logarithmic compression here
				// because we want to maintain visibility over wide dymanic range.
				res /= MAX_RES;
				res = 64 * log( 54.6 * res + 1 );
				get_colour( A->bits, &r, &g, &b );
				r *= res;
				g *= res;
				b *= res;
				r = floor( r + 0.5 );
				g = floor( g + 0.5 );
				b = floor( b + 0.5 );
				if( r < 0 ) r = 0;	if( r > 255 ) r = 255;
				if( g < 0 ) g = 0;	if( g > 255 ) g = 255;
				if( b < 0 ) b = 0;	if( b > 255 ) b = 255;
			}
			else
			{
				r = 0;
				g = 0;
				b = 0;
			}
			fputc(( int )r, fp );
			fputc(( int )g, fp );
			fputc(( int )b, fp );
			A++;
		}
	}
	fclose( fp );
	fprintf( stderr, "frame=%d\n", frame ); fflush( stdout ); fflush( stderr );
}

/**
 * Setup all the cell values to a reasonable starting value, with decent resources
 * and a flat tolerance function (i.e. ignore any similarity/difference) at a neutral
 * level (i.e. all players take a neutral attitude to each other).
 *
 * Then loop with random pairs (fast random generation is the bottleneck)
 * of neighbours, and transaction between the neighbours on the basis of above.
 * In addition, small random changes happen to the variables, so that mutation
 * exists in the system as a whole.
 */

int main()
{
	// Initial setup
	{
		int j, k;
		cell_T *A = grid;

		for( j = 0; j < HEIGHT; j++ )
		{
			for( k = 0; k < WIDTH; k++ )
			{
				setup( A );
				A++;
			}
		}
	}

	for( cycle_count = 0; cycle_count < TOTAL_CYCLES; cycle_count++ )
	{
		cell_T *A = 0;
		cell_T *B = 0;

		point_select( &A, &B );
		transaction( A, B );
		if(( cycle_count & 0x3FFFFF ) == 0 ) output_frame();
	}

}

#endif