// ecoNode.c
// K. J. Rock
// (c) August 13, 2024

// g++ -o en ecoNode.c
// gnuplot script
// set logscale y    or set logscale   for a log,log chart
// plot 'dt.dat' using 1:4 title 'Carnivore' with line, 'dt.dat' using 1:3 title 'Herbivore' with line, 'dt.dat' using 1:2 title 'Plant' with line
// plot 'dt.dat' using 2:3 with line, 'dt.dat' using 3:4 with line, 'dt.dat' using 2:4 with line
// plot 'dt.dat' using 1:4 title 'Carnivore' with line, 'dt.dat' using 1:3 title 'Herbivore' with line, 'dt.dat' using 1:2 title 'Plant' with line, 'dt.dat' using 1:5 with line, 'dt.dat' using 1:6 with line, 'dt.dat' using 1:7 with line
// Profile script:   g++ -o en -pg ecoNode.c  ./en creates gmon.out for gprof en

// https://www.quantamagazine.org/teen-mathematicians-tie-knots-through-a-mind-blowing-fractal-20241126/

// https://www.digitalocean.com/community/tutorials/how-to-enable-remote-desktop-protocol-using-xrdp-on-ubuntu-22-04

// curl ifconfig.me
// gives me my public ip address:: 208.67.206.202
// sudo ufw allow from 208.67.206.202/32 to any port 3389

#include <time.h>
#include <stdio.h>
#include <stdlib.h>

#define SZ 5000	// Height and Width are equal for a square 'toroidal' ecosystem
typedef unsigned int chromosome;					// Chromosome - 32 bits of genetic information
struct vt { int x, y; };							// 2D integer space
struct vt dir[9];										// 2D direction vector from D2Q9 LBM now mathematical
char G[32];												// Binary gene buffer for geneLog() & showGenes()

struct pnode		// Plant node
	{
	struct vt X;										// X is an i,j pair based on an SZxSZ grid
	float energy;										// Life force
	bool alive;											// Deva  
	int age;												// Metabolize
	chromosome EN, PL;								// EN for energy genes,  PL for plant genes
	chromosome GM;										// GM is a genetic marker passed directly from mom
	};

struct hnode		// Herbivore node
	{
	struct vt X;										// Location	
	int Psense[9], Csense[9];						// Sensor network
	float energy;										// Life force
	bool alive;											// Starve or get eaten?
	int age;												// Metabolize, entropy, senescence
	chromosome EN, HB;								// Energy and herbivore chromosomes
	chromosome SN, DR;								// Sensor and direction chromosomes
	};

struct cnode		// Carnivore node
	{
	struct vt X;										// 2D location vector
	int Hsense[9];										// Sensor network
	float energy;										// Chi
	bool alive;											// Energy above zero.
	int age;												// Used in reproduce()
	chromosome EN, CN;								// Energy and carnivore chromosomes
	chromosome SN, DR;								// Sensor and direction chromosomes
	};

struct enode		// Ecology node
	{
	float nutrient;									// The nutrient energy level of this cell
	struct pnode *P;									// Point to plant,  
	struct hnode *H;									// herbivore, and/or
	struct cnode *C;									// carnivore living here.
	};

	// Ecosystem and organism buffers  /////////
struct enode *E 	= (struct enode *) malloc( SZ * SZ * sizeof( struct enode ));
struct pnode *PL1 = (struct pnode *) malloc( SZ * SZ * sizeof( struct pnode ));
struct pnode *PL2 = (struct pnode *) malloc( SZ * SZ * sizeof( struct pnode ));
struct hnode *HL1 = (struct hnode *) malloc( SZ * SZ * sizeof( struct hnode ));
struct hnode *HL2 = (struct hnode *) malloc( SZ * SZ * sizeof( struct hnode ));
struct cnode *CL1 = (struct cnode *) malloc( SZ * SZ * sizeof( struct cnode ));
struct cnode *CL2 = (struct cnode *) malloc( SZ * SZ * sizeof( struct cnode ));

	// Lists for new plants, herbivores, and carnivores
struct pnode *NP = (struct pnode *) malloc( SZ * SZ * sizeof( struct pnode ));
struct hnode *NH = (struct hnode *) malloc( SZ * SZ * sizeof( struct hnode ));
struct cnode *NC = (struct cnode *) malloc( SZ * SZ * sizeof( struct cnode ));

	// Buffer indirection indices  /////////////
struct pnode *PL = PL1;								// Initial condition
struct hnode *HL = HL1;
struct cnode *CL = CL1;
struct pnode *OPL = PL2;							// Pointers to the "Other" buffers
struct hnode *OHL = HL2; 
struct cnode *OCL = CL2;
bool even = true;										// Buffer indirection flag

int newPpop, newHpop, newCpop;					// New plant, herbivore, carnivore counter
float pmin, pmax, hmin, hmax;						// Current high and low enery bounds
float cmin, cmax, emin, emax;						// Population extrema too
int Peat, Heat, Ceat;								// Number of eaten each cycle, or died

	// Main functions  /////////////////////////
void initEcosystem( void );
void initPlants( void );
void initHerbivores( void );
void initCarnivores( void );
void moveHerbivores( void );						// Call move() on each herbivore
void moveCarnivores( void );						// Call move() on each carnivore
float metabolize( void );							// metabolize()							
void replenish( float );							// replenish() nutrient layer in the E buffer
void eat( void );										// Absorb nutrients, eat plants, or eat herbivores
void reproduce( void );								// reproduce() if you have the energy and age
void die( void );										// die() from energy loss
void cull( void );									// Bring out your dead
void addNew( void );									// Catenate new organisms to their lists
void bufferSwap( void );							// Swap even/odd buffers

	// Useful tools  ///////////////////////////
struct vt vector( int, int );						// 2D integer vector for cell or pixel use
struct vt sum( struct vt, struct vt );			// Add 2D vectors
struct vt bound( struct vt );						// Place 2D vector into toroidal space
struct vt scale( int, struct vt );				// Find an open cell around current location
float fRan( void );									// Floating point random number generator
void initDir( void );								// Initialize the direction vector buffer
struct vt move( chromosome );						// Use motion genes to relocate animals
chromosome crossover( chromosome, chromosome );
chromosome mutation( chromosome );

	// Animal sensor array functions  //////////
void herbScan( void );								// Local sensor scan for plants and carnivores
void carnScan( void );								// Local sensor scan for herbivores
int maxHerb( int CS[] );							// Find best feeding vector

	// Data Logging functions  /////////////////
void storeData( void );								// Backup ecosystem to file
void loadData( void );								// Restore ecosystem from file
void storeSamples( void );							// Collect 1000 P,H,C samples
void mergeSamples( void );							// Add marked samples to eco.dat ecosystem
void loadSamples( void );							// Read the samples into their own buffers
void loadSampTwo( void );							// Read second sample into secondary buffers
void useSamples( void );							// Build an ecosystem from those samples
void showNew( void );								// Display data for last ten plants in PL
void showSamples( void );							// Display chromosomes of 1000 samples
void stamp( void );									// Add a time stamp output
void fStamp( FILE *fp );							// Stamp the filed data too

	//  Diagnostic tools  //////////////////////
void display( chromosome );						// Print chromosome binary view
void printPlants( void );							// Display a map of plants P or .
void printHerbivores( void );						// Display grid of herbivores H or .
void printCarnivores( void );						// Display grid of carnivores C or .
void checkOffspring( void );						// Display new offspring buffers
void check( void );									// Sanity checking code
void checkDir( void );								// Display direction vectors
void checkAlternate( void );						// Lucidity approaches
void checkEnergy( void );							// Super sanity assured
void checkOthers( void );							// Nirvana achieved
void geneLog( chromosome );						// Binary gene buffer
void showGenes( void );								// Show genes of 10 P, 10 H, 10 C
//void logEcosystem( void );

#ifdef NOTES
void logEcosystem( void )							// Initialize primary ecosystem buffer
	{
	FILE *fp = fopen("nut.dat", "w");			// Try storing less data
															// A 1000*1000 window inside of SZ*SZ?
	for(int j=0; j<SZ; j++)							// We could store palette index instead
		{
		for(int i=0; i<SZ; i++)		
			fprintf( fp, "%f", (E +j*SZ +i)->nutrient );
		fprintf( fp, "\n" );							// Terminate each line	
		}
	fclose( fp );
	}
#endif

int Ppop = 400000, Hpop = 200000, Cpop = 25000;	// Initial populations	

//#define FIRST										// First run of ecosystem
			// Necessary if you don't have: 'eco.dat', 'sample.dat', or 'sample9.dat'
int main( void )
	{
	int kk = 0;
	FILE *fp = fopen("dt.dat", "w");				// Instantiate data file
	fclose( fp );										// Clear data file for writing

	srand( (unsigned) time( NULL ) );			// Seed random number generator
	initDir();											// Fill dir[] array with motion vectors
	
#ifdef FIRST									// Create a new ecosystem with random values?
	initEcosystem();									// Establish ecosystem framework	
	initPlants();										// Create plants, fill ecosystem
	initHerbivores();									// Create herbivores, fill ecosystem
	initCarnivores();									// Create carnivores, fill ecosystem
#else												// Read ecosystem data from a file?
	loadData();											// Load ecosystem image file	'eco.dat'  100

//	loadSamples();										// Load the genetic sample file 'samples.dat'	
//	loadSampTwo();										// Load 'samples9.dat' for a second population
//	useSamples();										// Build a new ecosystem from samples  50
	
//	mergeSamples();									// Merge 'sample.dat' into 'eco.dat' system
//	showNew();											// Display data for last ten plants in PL 
#endif				

//	for(int zq = 0; zq < 20; zq++)				// Limit run for profiling with gprof
	while( true )										// Enter a forever loop
		{											// 'dt.dat' for population plotting data every cycle
//		printf("%4d %8d %7d %6d    %7d %6d %5d\n", kk, Ppop, Hpop, Cpop, Peat, Heat, Ceat );
		printf("%4d %8d %7d %6d    %6.2f %6.2f %6.2f   ", kk, Ppop, Hpop, Cpop, pmax, hmax, cmax );
		stamp();											// Add a time stamp to display
		fp = fopen("dt.dat", "a");					// Append to end of even a blank file
		fprintf(fp, "%d %d %d %d    %d %d %d\n", kk++, Ppop, Hpop, Cpop, Peat, Heat, Ceat );
//		fStamp( fp );									// Add time stamp to filed data
		fclose( fp );
//		logEcosystem();								// Store the nutrient layer for display.
		
		if( kk%250 == 0 )								// Store genetic samples
			storeSamples();							// in 'sample.dat' every 250 cycles
		if( kk%500 == 0 )								// Simulation image every 500 cycles
			storeData();								// in 'eco.dat' image of simulation
		
//		herbScan();										// Populate sensor grid  out 4  SLOW!!
		moveHerbivores();								// Check randomness
//		carnScan();										// Fill sensor array		 out 8
		moveCarnivores();								//	dir[] and scale()
		replenish( metabolize() );					// Balance energy out with energy in		
		eat();											// Surviving organisms gain energy
		reproduce();									// Mate organisms for the future
		die();											// Set organism alive flag to false :(
		cull();											// NULL the dead while collecting survivors
		addNew();										// Add offspring to 'other' buffers
		bufferSwap();									// Swap to 'other' buffers
		
		if( Ppop < 2 || Hpop < 2 || Cpop < 2 )	// It takes two to tango
			break;
		}
	}

//////////////////////////////////////////////////////////////////////////////////////////
//////  Functions()  /////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////

void initEcosystem( void )							// Initialize primary ecosystem buffer
	{
	for(int j=0; j<SZ; j++)
		for(int i=0; i<SZ; i++)
			{			
			(E +j*SZ +i)->nutrient = 250;			// PFA number
			(E +j*SZ +i)->P = NULL;					// No one lives here at present
			(E +j*SZ +i)->H = NULL;
			(E +j*SZ +i)->C = NULL;					// Dave's not here man
			}
	}

void initPlants( void )
	{
	int X, Y;											// Prospective location
	
	for(int i=0; i<Ppop; i++)						// Create a population of plants.
		{
		do	{						// Unless there is a living plant in that location.
			X = rand()%SZ;								// Generate new point.
			Y = rand()%SZ; 
			} 	while( (E +Y*SZ +X)->P );			// When a plant is resident
		
		(PL+i)->X.x = X;								// Remember where I live.	
		(PL+i)->X.y = Y;
		(E +Y*SZ +X)->P = (PL+i);					// Point at plant.
			
		(PL+i)->EN = (chromosome) rand();		// Create energy chromosome.
		(PL+i)->PL = (chromosome) rand();		// Create plant chromosome.
	
		if( fRan() < 0.07 )	
			(PL+i)->GM = 3855;						// Maternal genetic marker.
		else
			(PL+i)->GM = 0;
		(PL+i)->alive = true;						// I can feel the sun.
		(PL+i)->energy = 80;							// PFA number
		(PL+i)->age = 1;								// Day of birth		
		}
	}

void initHerbivores( void )
	{	
	int X, Y;											// Prospective location of herbivore
		
	for(int i=0; i<Hpop; i++)
		{	
		do	{												// Find an open location
			X = rand()%SZ;								// Generate prospective points
			Y = rand()%SZ; 
			} while( (E +Y*SZ +X)->H );			// Occupado
		
		(HL+i)->X.x = X;								// Where are my keys?
		(HL+i)->X.y = Y;	
		(E +Y*SZ +X)->H = (HL+i);					// Point at herbivore
			
		(HL+i)->EN = (chromosome) rand();		// Create energy chromosome
		(HL+i)->SN = (chromosome) rand();		// Create sensor chromosome
		(HL+i)->HB = (chromosome) rand();		// Create herbivore chromosome
		(HL+i)->DR = (chromosome) rand();		// Create direction chromosome
			// A marker chromosome would fit here  GM for genetic marker
		(HL+i)->alive = true;
		(HL+i)->energy = 120;						// PFA number
		(HL+i)->age = 1;
		}
	}

void initCarnivores( void )
	{
	int X, Y;											// Prospective location of herbivore
		
	for(int i=0; i<Cpop; i++)
		{				
		do	{
			X = rand()%SZ;								// Generate prospective points
			Y = rand()%SZ; 
			}  while( (E +Y*SZ +X)->C );			// Find your own lair
		
		(CL+i)->X.x = X;								// Remember where I live	
		(CL+i)->X.y = Y;
		(E +Y*SZ +X)->C = (CL+i);  				// Point at carnivore

		(CL+i)->EN = (chromosome) rand();		// Create energy chromosome
		(CL+i)->SN = (chromosome) rand();		// Create sensor chromosome
		(CL+i)->CN = (chromosome) rand();		// Create carnivore chromosome
		(CL+i)->DR = (chromosome) rand();		// Create direction chromosome
			// A marker chromosome would fit here  GM for genetic marker
		(CL+i)->alive = true;
		(CL+i)->energy = 300;						// PFA number
		(CL+i)->age = 1;
		}
	}
	
void moveHerbivores( void )
	{
	struct vt A, B, C;								// Working vectors
	int min, ring, ch, speed;

	for(int i=0; i<Hpop; i++)						// Scan herbivore list
		{
		ch = 0;											// Limit choice count
		min = 1;											// { min <= speed <= 8 }		
		speed = (int) (((HL+i)->HB & (7 << 28)) >> 28 ) + min;
		
		A = (HL+i)->X;									// Current location		
		do	{												// Find an open spot within range	
			ring = min + rand()%speed;				// {min <= ring < speed}		
			B = scale( ring, move( (HL+i)->DR ) );	// Express direction genes		
			C = bound( sum( A, B ) );				// Prospective new location			
			if( ch > 8 )  break;  else  ch++;	// Limit search count to 8 choices
			} 	while( (E +C.y*SZ +C.x)->H );		// While location is occupied
		
		(E +A.y*SZ +A.x)->H = NULL;				// Empty old location
		(HL+i)->X = C;									// Remember where I live		
		(E +C.y*SZ +C.x)->H = (HL+i);				// Fill new location			
		}
	}

void herbScan( void )
	{
	struct vt A, B, C;								// Working vectors
	int maxRing = 3;  								// PFA number
	
	for(int i=0; i<Hpop; i++)
		{
		A = (HL+i)->X;									// Current location		
		for(int k=1; k<9; k++)						// Clear Psense[] and Csense
			{
			(HL+i)->Psense[ k ] = 0;
//			(HL+i)->Csense[ k ] = 0;
			}
			// Scan from inner ring to outer ring defined by sensor strength
		for(int ring=1; ring<maxRing; ring++)			
			for(int k=1; k<9; k++)					// Scan directions
				{				
				B = scale( ring, dir[ k ] );  	// Walk around ring 1..8
				C = bound( sum( A, B ) );			// Search location
				if( (E +C.y*SZ +C.x)->P)  (HL+i)->Psense[ k ]++;	// I found a plant
//				if( (E +C.y*SZ +C.x)->C)  (HL+i)->Csense[ k ]++;	// I found a carnivore
				}
		}
	}

void carnScan( void )								// Sense herbivores in local area
	{														// 8*8 = 64 possible cells
	int maxRing = 8;									// PFA number	potential gene  111 + 1								
	struct vt A, B, C;								// Working vectors
	
	for(int i=0; i<Cpop; i++)
		{
		A = (CL+i)->X;									// You are here.
										
		for(int k=1; k<9; k++)						// Clear Hsense[]
			(CL+i)->Hsense[ k ] = 0;
			
		for(int ring=1; ring<maxRing; ring++)
			for(int k=1; k<9; k++)
				{
				B = scale( ring, dir[ k ] );		// Walk out the direction vector
				C = bound( sum( A, B ) );			// Search location				
				if( (E +C.y*SZ +C.x)->H )  (CL+i)->Hsense[ k ]++;	// I found a herbivore :)
				}
		}
	}

#ifdef NOTES
Carnivores know their herbivore environment through their individual sensor arrays.
Find the largest number in the array[ 1..8 ]
Then scale/step by ring along that vector, feasting on herbivores.
Continue until array[ k ] count is done, or you have run out of hunting choices
Why not decrement array[ k ]-- as Herbivores are eaten under its direction.
Then when array[ k ] == 0 you can increment k
Or find the the remaining highest count.
Repeat stuff until the choice count is exceeded.
Try to make this quicker than the herbScan() and carnScan() functions.

CS[ i ] gives us direction and quantity
If we decrement the quantity after eating a herbivore
We can call maxHerb() again for the next best feeding vector
repeat until XX are eaten  OR we need another break out

Where should this be implemented?  In the eat() function?
Or should it be in the move() function?
Or in moveCarnivores()??  I am confused :)

In eat() I could sharpen my focus and eat along high yield vectors.
If I put it into moveCarnivores() I could direct them to where their 
sensor arrays aim them.					TODO

We could make parameters of the sensor array genetic
Or make it cost more energy to use herbScan() over random choices.   ???

speed = (int) (((CL+i)->CN & (15 << 27)) >> 27) + 12;	// CN::G1	12 - 27
ring = rand()%speed;						// {min <= ring <= speed}	DR::G1..G8
B = scale( ring, move( (CL+i)->DR ) );	// Motion vector from DR chromosome	
C = bound( sum( A, B ) );				// Prospective new location

For 57 random choices
#endif

// Script to back up the working BD repository:
// svnrdump dump svn+ssh://hermes/srv/samba/share3/BD_2nd_Edition | gzip -9 > ~/TEST/repo.dump.gz

int maxHerb( int CS[] )								// 
	{
	int Hmax;											// Most dense direction
	int max = -100;									// Guard
	for(int i=1; i<9; i++)
		if( CS[ i ] > max )							// Compare with current maximum
			{
			Hmax = i;									// Seek supremum
			max = CS[ Hmax ]; 						// One shot 
			}												// Recaluculate each step ??
	return Hmax;
	}

void moveCarnivores( void )
	{
	struct vt A, B, C;								// Current position, move, prospective
	int min, ring, ch, speed;

	for(int i=0; i<Cpop; i++)						// Scan carnivore list
		{
		ch = 0;											// Choice counter
		min = 3;											// { min <= speed <= 10 }
		speed = (int) (((CL+i)->CN & (7 << 28)) >> 28) + min;
				
		A = (CL+i)->X;									// Current location				
		do	{												// Find an open spot within range
			ring = min + rand()%speed;				// {min <= ring < speed}
			B = scale( ring, move( (CL+i)->DR ) );	// Motion vector from DR chromosome		
			C = bound( sum( A, B ) );				// Search location	

			if( ch > 8 ) break;  else  ch++;		// All your choice are us
			}  while( (E +C.y*SZ +C.x)->C );		// Carnivore already lives here

		(E +A.y*SZ +A.x)->C = NULL;				// Blank old spot
		(CL+i)->X = C;									// Remember where I live
		(E +C.y*SZ +C.x)->C = (CL+i);				// Inhabit new spot
		}
	}

////////////////////////////////////////////////////////////////////////////////////////
/////////////  Parameter Tweaking Section  /////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////// 
	
float metabolize( void )  							// Sum energy lost by all organisms
	{
	float energy = 0;  								// Balance global energy loss and gain 
	
	for(int i=0; i<Ppop; i++)  					// EN::G5 energy use per cycle
		{				
		(PL+i)->age += 1;								//	EN::G5  7.5 - 14.5
		(PL+i)->energy -= (float) ((PL+i)->EN & 7) + 7.5;	// PFA number // TODO 6.1
		energy += (float) ((PL+i)->EN & 7) + 7.5;				// PFA number
		}
	for(int i=0; i<Hpop; i++)
		{
		(HL+i)->age += 1;								//	EN::G5  4.2 - 11.2
		(HL+i)->energy -= (float) ((HL+i)->EN & 7) + 4.2;	// PFA number // TODO 4.2
		energy += (float) ((HL+i)->EN & 7) + 4.2;				// PFA number
		}
	for(int i=0; i<Cpop; i++)
		{
		(CL+i)->age += 1;			 					//	EN::G5  5.3 - 20.3
		(CL+i)->energy -= (float) ((CL+i)->EN & 7) + 5.3;	// PFA number TODO 5.3
		energy += (float) ((CL+i)->EN & 7) + 5.3;
		}
	return energy;										// Pass energy use to replenish()
	}

void replenish( float energy )					// Balance energy lost during metabolism
	{
	emin = 1e5;		emax = -1e5;					// Guard the endpoints
	
	for(int i=0; i<SZ*SZ; i++)
		{													// PFA number TODO 2.7 1.5
		(E+i)->nutrient += energy / (SZ*SZ) + 2.3;	// Spread energy evenly across ecosystem
															// TODO kjr November 7th, 2024
		if( (E+i)->nutrient < emin ) emin = (E+i)->nutrient;			// Determine extrema
		if( (E+i)->nutrient > emax ) emax = (E+i)->nutrient;
		}
	}

void die( void )
	{
	int senescence;									// EN::G1	
	Ceat = 0;											// Clear carnivore eaten counter							
	
	for(int i=0; i<Ppop; i++)
		{													// 2 - 5
//		senescence = (int) (((PL+i)->EN & (3 << 22)) >> 22) + 2;		// PFA number
		senescence = 5;								// PFA number
		if( ((PL+i)->energy < 0) || ((PL+i)->age > senescence) )
			(PL+i)->alive = false;					// Like a tumblin' tumbleweed.
		}
	for(int i=0; i<Hpop; i++)
		{													// 19 - 50
		senescence = (int) (((HL+i)->EN & (31 << 22)) >> 22) + 19;	// PFA number
		if( ((HL+i)->energy < 0) || ((HL+i)->age > senescence) )
			(HL+i)->alive = false;					// Dave's not here man!
		}
	for(int i=0; i<Cpop; i++)
		{													// 39 - 70
		senescence = (int) (((CL+i)->EN & (31 << 22)) >> 22) + 39;	// PFA number
		if( ((CL+i)->energy < 0) || ((CL+i)->age > senescence) )
			{
			(CL+i)->alive = false;					// This parrot has ceased to be.
			Ceat++;		// Count Ceat here even though it's not really eaten, just dead.
			}
		}	
	}
															// Express thirty genes ??
void eat( void )										// Express fourteen genes
	{							
	struct pnode *P;									// Point at plants
	struct hnode *H;									// Point at herbivores
	struct vt A, B, C;								// Working vectors
	int j, ring;										// Direction, ring#
	float satiety, absorb;							// Genetic expression buffers	
	int ch, speed;										// Choice counter, express gene HB::G1
	float spine, eat;									// Plant vs Herbivore spine war
	float strength, withstand;						// Plant vs Herbivore poison war
	float camo, see;									// Herbivore vs Carnivore camo war
	
	Peat = Heat = 0;									// Count the eaten
															// Ceat is calculated in die()
		// Plants eat nutrients at E(x,y) only
	for(int i=0; i<Ppop; i++)						// Express two genes
		{
		C = (PL+i)->X;													
															// EN::G3  19.1 - 82.1
		satiety = (float) (((PL+i)->EN & (63 << 10)) >> 10) + 19.1;		// PFA number
		if( (PL+i)->energy < satiety )			// Am I hungry?
			{
//			printf("p");								// EN::G4  6.25 - 50.0%  1/16 to 8/16  
			absorb = (float) ((((PL+i)->EN & (7 << 7)) >> 7) + 1.0) / 15.0;
			(PL+i)->energy += (E +C.y*SZ +C.x)->nutrient * absorb;				// Add energy to plant			
			(E +C.y*SZ +C.x)->nutrient -= (E +C.y*SZ +C.x)->nutrient * absorb;// Remove from nutrient layer
			}			
		}	

		// Herbivores eat at E(x,y) and the rings around it	
	for(int i=0; i<Hpop; i++)						// Express fifteen genes
		{													// Spine wars PL::G2 vs HB::G4  PL::G3 vs HB::G3
		ch = 0;											// Limit choice count 			
		speed = (int) (((HL+i)->HB & (15 << 27)) >> 27) + 12;	// HB::G1   12 - 27 PFA number		
		eat = (float) (((HL+i)->HB & (31 << 15)) >> 15) / 31.0;	// HB::G4  ability to eat spines
		withstand = (float) (((HL+i)->HB & (31 << 15)) >> 15) / 31.0; // HB::G3 withstand poison
		satiety = (float) (((HL+i)->EN & (63 << 10)) >> 10) + 56.1;	// EN::G3 56.1 - 119.1

		A = (HL+i)->X;									// Current location					
		while( true )									// Forever loop requires break out
			{				// TODO should this be eating radius instead?  Express another Gene!!!		
			ring = rand()%speed;						// {min <= ring <= speed}	DR::G1..G8
			B = scale( ring, move( (HL+i)->DR ) );	// Express direction genes
			C = bound( sum( A, B ) );				// Prospective new location 
					
			if( (E +C.y*SZ +C.x)->P && ((HL+i)->energy < satiety ))
				{											// PFA number   63+31 = 94  63/94 = 0.67
				P = (E +C.y*SZ +C.x)->P;  				// Point to plant at E(x,y)  EN::G4  0.1 - 0.8 
				spine = (float) ((P->PL & (63 << 18)) >> 18) / 94.0;  // PL::G2  spine probability
				strength = (float) ((P->PL & (63 << 12)) >> 12) / 94.0; // PL::G3  poison strength
				if( (eat >= spine) && (withstand >= strength ))	
					{ 							// HL::HB::G4 > PL::PL::G2   &&   HL::HB::G3 >= PL::PL::G3					
					absorb = (float) ((((HL+i)->EN & (7 << 7)) >> 7) + 1.0) / 10.0; 	// PFA number
					(HL+i)->energy += P->energy * absorb;		// Add energy to herbivore			
					P->energy -= P->energy * absorb;				// Remove energy from plant
					(E +C.y*SZ +C.x)->nutrient += P->energy;	// Remainder to nutrient layer
					(E +C.y*SZ +C.x)->P = NULL;	// Plant no longer exists here
					P->alive = false;					// Plant is deceased	
					Peat++;								// Count the fodder	
					}
				}
							// Ring 12 to 27  8*12 = 96  27*8 = 216  So 32/96 or 32/216 is OK
			if( ch > 29 )  break;  else ch++;	// Limit search count		
			}			
		}
	
		// Carnivores eat in rings around E(x,y) along the 8 primary directions.
	for(int i=0; i<Cpop; i++)						// Express thirteen genes
		{		
		ch = 0;											// Choice counter
		speed = (int) (((CL+i)->CN & (15 << 27)) >> 27) + 12;	// CN::G1	12 - 27 
		see = (float) ((CL+i)->SN & 255) / 255.0;					// SN::G5  sight ability
		satiety = (float) (((CL+i)->EN & (7 << 10)) >> 10) + 8.0; // 8.0 - 15.0
															// TODO kjr November 7th, 2024
		A = (CL+i)->X;									// Remember location.		
		
		while( true )									// Eat until break
			{
			ring = rand()%speed;						// {min <= ring <= speed}	DR::G1..G8
			B = scale( ring, move( (CL+i)->DR ) );	// Motion vector from DR chromosome	
			C = bound( sum( A, B ) );				// Prospective new location
			
			if( (E +C.y*SZ +C.x)->H && ((CL+i)->energy < satiety) )
				{											// HB::G5  camouflage effectiveness
				H = (E +C.y*SZ +C.x)->H;			// Point to herbivore at E(x,y)	1/9 to 8/9
				camo = (float) ((H->HB & (31 << 10)) >> 10) / 50.0;   // Herbivore camouflage
				if( see >= camo )						// CL::SN::G5 >= HL::HB::G5
					{										// TODO kjr November 7th, 2024
					absorb = (float) ((((CL+i)->EN & (7 << 7)) >> 7) + 1.0) / 8.0;	// 0.125 - 1.0
					(CL+i)->energy += H->energy * absorb;	// Add energy to carnivore
					H->energy -= H->energy * absorb;			// Remove energy from herbivore				
					(E +C.y*SZ +C.x)->nutrient += H->energy;	// Remainder to nutrient layer					
					(E +C.y*SZ +C.x)->H = NULL;	// Herbivore no longer exists here					
					H->alive = false;					// Herbivore is deceased
					Heat++;								// Count the fallen							
					}
				}
						// Is 57 excessive?
			if( ch > 57 )  break;  else ch++;	// All your choice are us.
			}
		}	
	}														// End of eat()			

void reproduce( void )
	{
	int maturity, nw;									// Age, new organism counter
	float satiety;										// Am I hungry?
	struct vt A, B, C;								// Working vectors
	
	newPpop = newHpop = newCpop = 0;				// Clear population counters here
		// 4000 * 4000 = 16,000,000   9/16 = 0.5625  9.5 Million seg faulted
	if( Ppop < SZ*SZ*0.5 )							// HARD limit proportional to grid size
		{
		nw = 0;
		for(int i=0; i<Ppop; i++)
			{			
//			maturity = (int) (PL+i)->EN & 7 + 3;  	// EN::G4, ignore 2 high order bits	
//			satiety = (float) (((PL+i)->EN & (7 << 10)) >> 10) + 1.0;	EN::G3  1 - 512
//			if( ( (PL+i)->age > maturity ) && ((PL+i)->energy > satiety ) )	// TODO
			if( ((PL+i)->age > 3) && ((PL+i)->energy > 47.5) )			// PFA number 25.1
				{											// && pollen > 0.25???
				A = (PL+i)->X;							// Location of parent plant.			
				int r = rand()%Ppop;					// Select random mate 

				if( ((PL+r)->age > 3) && ((PL+r)->energy > 47.5) )		// PFA number
					{			
						// Need to find an empty spot in E
					int j = 1;							// Skip center spot
					int ring = 5;						// Pollination inner boundary
									
					do	{
						B = scale( ring, dir[j] );	// Motion vector
						C = bound( sum( A, B ) );	// Prospective new location	
						
						j++;								// Walk around direction vector					
						if( j > 8 ) 					// End of dir[]
							{
							j = 1;						// Search next
							ring++;						// outer ring
							}
						if( ring > 20 )  break;		// The Outer Limit!!
						} 	while( (E +C.y*SZ +C.x)->P );	// Occupied
						
						// New plant needs to be filled
					(NP+nw)->X = C;					// Put plant at new location				
					(NP+nw)->energy = (PL+i)->energy * 0.125;		// Natal energy
					(NP+nw)->energy += (PL+r)->energy * 0.125;	// 1/8 from each parent
			
					(PL+i)->energy -= (PL+i)->energy * 0.125;  	// Parent loses 1/8 of its energy
					(PL+r)->energy -= (PL+r)->energy * 0.125;		// Other parent loses energy
					
					(NP+nw)->EN = crossover( (PL+i)->EN, (PL+r)->EN );	// Crossover both plant's
					(NP+nw)->PL = crossover( (PL+i)->PL, (PL+r)->PL );	// Chromosomes
					(NP+nw)->GM = (PL+i)->GM;  	// Pass Mom's genetic marker on intact.
/*					if( (NP+nw)->GM == 3855 )		// Mutate marked plants
						if( fRan() > 0.5 )
							{
							(NP+nw)->EN = mutation( (NP+nw)->PL );
							(NP+nw)->PL = mutation( (NP+nw)->PL );
							}
*/						
															// My name is 905
					(NP+nw)->alive = true;        // and I've just become alive
					(NP+nw)->age = 1;					// I'm the newest populator  					
															// of the planet called Earth.
					(E +C.y*SZ +C.x)->P = (NP+nw);	// Point to plant data
					nw++;  								// Add new plant to count	
					}	
				}				
			}
		newPpop = nw;									// Set new plant list counter	
		}													// End of plant population limiter
	
	nw = 0;
	for(int i=0; i<Hpop; i++)
		{
		maturity = (int) (HL+i)->EN & 3 + 1;  	// EN::G4, ignore high order bit		
//		satiety = (float) (((HL+i)->EN & (31 << 10)) >> 10) + 21.0;
//		if( ( (HL+i)->age > maturity ) && ((HL+i)->energy > satiety ) )	// TODO				
		if( ((HL+i)->age > maturity) && ((HL+i)->energy > 105.0) )  // PFA number
			{
			A = (HL+i)->X;								// Location of parent plant.			
			int r = rand()%Hpop;						// Select random mate 
			
				// Check whether (PL+r) is Old enough
			maturity = (int) (HL+r)->EN & 3 + 1;	// At least one day old	
//			satiety = (float) (((HL+i)->EN & (31 << 10)) >> 10) + 21.0;
//			if( ( (HL+r)->age > maturity ) && ((HL+r)->energy > satiety ) )	// TODO
			if( ((HL+r)->age > maturity) && ((HL+r)->energy > 105.0) )  // PFA number 125.0
				{	// Need to find an empty spot in E
				int j = 1;								// Skip center spot
				int ring = 4;
					
				do	{
					B = dir[ j ];  					// Walk around each direction
					B = scale( ring, B );			// Motion vector
					C = bound( sum( A, B ) );		// Prospective new location	
					
					j++;									// Walk around direction vector					
					if( j > 8 )							// End of dir[]
						{
						j = 1;							// Search
						ring++;							// next outer ring
						}
					if( ring > 14 )  break;
					} 	while( (E +C.y*SZ +C.x)->H );	// Is this seat taken?

					// New herbivore needs to be prepared
				(NH+nw)->X = C;						// Remember where I left my keys
				(NH+nw)->energy = (HL+i)->energy * 0.125;		// Offspring gains natal energy
				(NH+nw)->energy += (HL+r)->energy * 0.125;	// 1/8 from each parent
		
				(HL+i)->energy -= (HL+i)->energy * 0.250;  	// Mother loses 1/4 of its energy
				(HL+r)->energy -= (HL+r)->energy * 0.125;		// Father loses 1/8 of its energy
				(E +C.y*SZ +C.x)->nutrient += (HL+i)->energy * 0.125;  // Afterbirth from mother				
				
				(NH+nw)->EN = crossover( (HL+i)->EN, (HL+r)->EN );	// Crossover: Energy
				(NH+nw)->DR = crossover( (HL+i)->DR, (HL+r)->DR );	// Direction genes
				(NH+nw)->SN = crossover( (HL+i)->SN, (HL+r)->SN );	// Sensor genes
				(NH+nw)->HB = crossover( (HL+i)->HB, (HL+r)->HB );	// Herbivore genes	
				// Genetic marker goes here									
				(NH+nw)->alive = true;				// The Quickening
				(NH+nw)->age = 1;  					// Day one of a new life							
				(E +C.y*SZ +C.x)->H = (NH + nw);	// Take up residence in my new home
				nw++;  									// Add new plant to count	
				}	
			}				
		}
	newHpop = nw;										// Set new herbivore list counter
	
	nw = 0;
	for(int i=0; i<Cpop; i++)						// Reproduce carnivores	
		{
		maturity = (int) (CL+i)->EN & 7 + 3;  	// EN::G4		
//		satiety = (float) (((CL+i)->EN & (15 << 10)) >> 10) + 20.0;
//		if( ( (CL+i)->age > maturity ) && ((CL+i)->energy > satiety ) )
		if( ((CL+i)->age > maturity) && ((CL+i)->energy > 39.5) )	// PFA number 43.7
			{												// TODO kjr November 7th, 2024
			A = (CL+i)->X;								// Location of parent herbivore.			
			int r = rand()%Cpop;						// Select random mate 
			
				// Check whether (PL+r) is Old enough
			maturity = (int) (CL+r)->EN & 7 + 3;  		// EN::G4, ignore high order bit		
//			satiety = (float) (((CL+i)->EN & (15 << 10)) >> 10) + 20.0;
//			if( ((CL+r)->age > maturity) && ((CL+r)->energy > satiety) )
			if( ((CL+r)->age > maturity) && ((CL+r)->energy > 39.5) )	// PFA number	
				{											// TODO kjr November 7th, 2024
					// Need to find a empty spot in E and in OE
				int j = 1;								// Skip the center spot
				int ring = 3;							// Inner loop
				
				do	{		
					B = dir[ j ];						// Walk around
					B = scale( ring, B );			// the motion vector
					C = bound( sum( A, B ) );		// Prospective new location
					
					j++;									// Walk around motion vector
					if( j > 8 )				  			// Limit search count
						{
						j = 1;
						ring++;
						}
					if( ring > 9 )  break;
					} 	while( (E +C.y*SZ +C.x)->C );

					// New carnivore needs to be filled
				(NC+nw)->X = C;						// Put carnivore at new location				
				(E +C.y*SZ +C.x)->C = (NC + nw);	// Point to plant data
				
				(NC+nw)->energy = (CL+i)->energy * 0.125;		// Natal energy
				(NC+nw)->energy += (CL+r)->energy * 0.125;	// 1/8 from each parent		
				(CL+i)->energy -= (CL+i)->energy * 0.25;  	// Mama loses 1/4 of its energy
				(CL+r)->energy -= (CL+r)->energy * 0.125;		// Papa loses 1/8 of its energy
				(E +C.y*SZ +C.x)->nutrient += (CL+i)->energy * 0.125;  // Afterbirth from mama
				
				(NC+nw)->EN = crossover( (CL+i)->EN, (CL+r)->EN );	// Crossover 
				(NC+nw)->SN = crossover( (CL+i)->SN, (CL+r)->SN );	// Chromosomes
				(NC+nw)->CN = crossover( (CL+i)->CN, (CL+r)->CN );	// for offspring
				(NC+nw)->DR = crossover( (CL+i)->DR, (CL+r)->DR );	// genetics
				// Genetic marker goes here
				(NC+nw)->alive = true;
				(NC+nw)->age = 1;  					// Newly born				
						
				nw++;  									// Add new plant to count	
				}	
			}			
		}
	newCpop = nw;										// Set new plant list counter	with offset
	}														// End of reproduce()	

void cull( void )										// Garbage collection using twin buffers
	{
	int k=0, q=0, c=0;								// Organism counters  Need globally
	struct vt T;										// Working vector	
	pmin = hmin = cmin = 1e5;						// Guard the Outer Limits
	pmax = hmax = cmax = -1e5;
	
	for(int i=0; i<Ppop; i++)
		if( (PL+i)->alive )
			{
			if( (PL+i)->energy < pmin ) pmin = (PL+i)->energy;
			if( (PL+i)->energy > pmax ) pmax = (PL+i)->energy;
			*(OPL+k) = *(PL+i);  					// Point at survivor node from alternate buffer										
			k++;											// Count the survivors.
			}
		else
			{
			T = (PL+i)->X;								// Index Ecosystem plant location
			(E +T.y*SZ +T.x)->P = NULL;			// NULL plant pointer 
			}			
					
	for(int i=0; i<Hpop; i++)
		if( (HL+i)->alive )
			{
			if( (HL+i)->energy < hmin ) hmin = (HL+i)->energy;
			if( (HL+i)->energy > hmax ) hmax = (HL+i)->energy;
			*(OHL+q) = *(HL+i);						// Does indeed deep copy all information										
			q++;											// Not quite, it is simply that the pointers
			}												// Point to the same data.  
		else
			{			
			T = (HL+i)->X;
			(E +T.y*SZ +T.x)->H = NULL;	
			}
			
	for(int i=0; i<Cpop; i++)						// There is no copying involved.
		if( (CL+i)->alive )
			{
			if( (CL+i)->energy < cmin ) cmin = (CL+i)->energy;
			if( (CL+i)->energy > cmax ) cmax = (CL+i)->energy;
			*(OCL+c) = *(CL+i);						// Copy survivor to 'other' list							
			c++;											// Count survivors
			}
		else
			{		
			T = (CL+i)->X;
			(E +T.y*SZ +T.x)->C = NULL;			
			}
 
	Ppop = k;  Hpop = q;  Cpop = c;				// Update counters
	}														// End of cull()

void addNew( void )
	{
	int X, Y;

	for(int i=0; i<newPpop; i++)
		{
		*(OPL+Ppop+i) = *(NP+i);					// Copy plant into 'other' buffer
		X = (NP+i)->X.x;								// Create indices
		Y = (NP+i)->X.y;	
		(E +SZ*Y +X)->P = (OPL+Ppop+i);			// Add new plant to ecosystem
		}
	
	for(int i=0; i<newHpop; i++)
		{
		*(OHL+Hpop+i) = *(NH+i);					// Copy herbivore into 'other' buffer
		X = (NH+i)->X.x;
		Y = (NH+i)->X.y;
		(E +SZ*Y +X)->H = (OHL+Hpop+i);			// Add new herbivore to ecosystem
		}
	
	for(int i=0; i<newCpop; i++)
		{
		*(OCL+Cpop+i) = *(NC+i);					// Copy into 'other' buffer
		X = (NC+i)->X.x;
		Y = (NC+i)->X.y;
		(E +SZ*Y +X)->C = (OCL+Cpop+i);			// Add new carnivore to ecosystem
		}
	
	Ppop += newPpop;  Hpop += newHpop;  Cpop += newCpop;
	}														// End of addNew()

void bufferSwap( void )
	{
	if( even )
		{  
		PL  = PL1; 	HL = HL1;   CL = CL1;		// Even step buffers
		OPL = PL2; OHL = HL2;  OCL = CL2;		// 'Other' list pointers :)
		even = false;
		}
	else
		{
		PL  = PL2; 	HL = HL2;  	CL = CL2;		// Odd step buffers
		OPL = PL1; OHL = HL1;  OCL = CL1;
		even = true;
		}
	}

////////////////////////////////////////////////////////////////////////////////////////
/////////////  Useful Tools  ///////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////// 

float fRan( void )
	{
	return (float) rand() / (float) RAND_MAX;
	}

struct vt vector( int x, int y )
	{
	struct vt A;
	A.x = x;  A.y = y;
	return A;
	}

struct vt sum( struct vt A, struct vt B )
	{
	struct vt C;
	C.x = A.x + B.x;
	C.y = A.y + B.y;
	return C;
	}		

struct vt scale( int sc, struct vt B )
	{
	struct vt A;
	A.x = B.x * sc;
	A.y = B.y * sc;
	return A;
	}

struct vt bound( struct vt C )
	{
	if( C.x < 0 )				C.x += SZ;	// Connect ends
	else if( C.x > SZ-1 ) 	C.x -= SZ;	// to wrap 
	if( C.y < 0 )  			C.y += SZ;	// around into
	else if( C.y > SZ-1 )  	C.y -= SZ;	// a toroidal space.
	return C;
	}

	// Mathematical dir[]  :: { E, NE, N, NW, W, SW, S, SE }
void initDir( void )
	{	
	int x = 1, y = 1;									// For unit vectors
	dir[0] = vector( 0, 0 );						// Center
	dir[1] = vector( x, 0 );						// East
	dir[2] = vector( x, y );						// Northeast
	dir[3] = vector( 0, y );						// North
	dir[4] = vector( -x, y );						// Northwest
	dir[5] = vector( -x, 0 );						// West
	dir[6] = vector( -x, -y );						// Southwest
	dir[7] = vector( 0, -y );						// South
	dir[8] = vector( x, -y );						// Southeast
	}

struct vt move( chromosome DR )  				// Get direction from chromosome
	{	
   float dirProb[9];  								//		LBM version			Math version
	dirProb[8] = (float)  (DR & 15);							// SW		DR::G8	SE
	dirProb[7] = (float) ((DR & (15 << 4)) >> 4);		// NW					S
	dirProb[6] = (float) ((DR & (15 << 8)) >> 8);		// SE					SW
	dirProb[5] = (float) ((DR & (15 << 12)) >> 12);		// NE					W
	dirProb[4] = (float) ((DR & (15 << 16)) >> 16);		// S					NW
	dirProb[3] = (float) ((DR & (15 << 20)) >> 20);		// N					N
	dirProb[2] = (float) ((DR & (15 << 24)) >> 24);		// W					NE
	dirProb[1] = (float) ((DR & (15 << 28)) >> 28);		// E		DR::G1	E
	dirProb[0] = 0.0;									// 0 is at the center 
	
	float tt = 0;  									// Eight directions, 4 bit probability range
	for(int i=1; i<9; i++)
		tt += dirProb[i];								// Sum probabilities
	
	for(int i=1; i<9; i++)
		dirProb[i] /= tt;								// Normalize direction probabilities

	float rn = fRan();  								// Gives us a number from 0.0 to 1.0 inclusive
	float sum = 0;										// Sum probabilities until they exceed rn 
	int dr = 0;											// Direction buffer
   while (sum < rn)
   	{
      dr++;												// Sum probabilities from 1 to 8 if necessary
      if (dr > 8) dr = 0;
      sum += dirProb[ dr ];						// Probabilities must total 1.00
      }													// for this to work properly

	return dir[ dr ];									// Prospective direction
   }

chromosome crossover( chromosome CM, chromosome CP ) // Chromosome crossover
	{
	int i = rand() % 31;								// Determine crossover point
	
	CM >>= i;								      	// Clear upper bits of Mama chromosome
	CM <<= i;
	CP <<= 31 - i;						         	// Clear lower bits of Papa chromosome
	CP >>= 31 - i;
	
	if( fRan() > 0.95 )
		return mutation( CM | CP );				// Mutate chromosome 5.0% of the time
	else
		return ( CM | CP );			            // Merge upper and lower parts
	}

chromosome mutation( chromosome C )		   	// Mutate DNA in this chromosome
	{
	int i = rand() % 31;			
	return C ^= (1 << i);							// Flip a single bit at i
	}

/////////////////////////////////////////////////////////////////////////////////////////////////
///////////////  Data Logging Functions  ////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////

void storeData( void )
	{
	FILE *fp;
//	fp = fopen("eco.dat", "w");	
	static int q = 0;		// Snapshot number
	char buf[10];
	sprintf( buf, "eco%02d.dat", q++ );
	fp = fopen( buf, "w" );
															// Store SZ, Ppop, Hpop, Cpop;
	fprintf(fp, "%d %d %d %d\n", SZ, Ppop, Hpop, Cpop );		
	
	for(int i=0; i<SZ*SZ; i++)						// Store nutrient layer
		fprintf(fp, "%f ", (E+i)->nutrient );
	fprintf(fp, "\n");
	
	for(int i=0; i<Ppop; i++)						// Store only living organisms :) 
		fprintf(fp, "%d %d %u %u %u %f %d\n", 	(PL+i)->X.x, (PL+i)->X.y, (PL+i)->EN, (PL+i)->PL, 
														(PL+i)->GM, (PL+i)->energy, (PL+i)->age );	
	for(int i=0; i<Hpop; i++)
		fprintf(fp, "%d %d %u %u %u %u %f %d\n", (HL+i)->X.x, (HL+i)->X.y, (HL+i)->EN, (HL+i)->SN,
														(HL+i)->HB, (HL+i)->DR, (HL+i)->energy, (HL+i)->age );
	for(int i=0; i<Cpop; i++) 
		fprintf(fp, "%d %d %u %u %u %u %f %d\n", (CL+i)->X.x, (CL+i)->X.y, (CL+i)->EN, (CL+i)->SN,
														(CL+i)->CN, (CL+i)->DR, (CL+i)->energy, (CL+i)->age );	
	fclose( fp );
	}

void loadData( void )
	{
	int SX;  											// Temporary for SZ
	FILE *fp;
	fp = fopen("eco.dat", "r");
//	fp = fopen("eco01.dat", "r");	
															// Read parameters		
	fscanf( fp, "%d %d %d %d", &SX, &Ppop, &Hpop, &Cpop );	
	
	for(int i=0; i<SX*SX; i++)						// Refill the nutrient layer
		fscanf(fp, "%f", &(E+i)->nutrient );
	
	for(int i=0; i<Ppop; i++)						// Refill the plant buffer
		{
		fscanf(fp, "%d %d %u %u %u %f %d", &(PL+i)->X.x, &(PL+i)->X.y, &(PL+i)->EN, &(PL+i)->PL, 
											 		  &(PL+i)->GM, &(PL+i)->energy, &(PL+i)->age );		
/*		if( fRan() < 0.27 )							// 27% marking rate
			(PL+i)->GM = 3855;						// Maternal genetic marker
		else
			(PL+i)->GM = 0;		*/
		}
	
	for(int i=0; i<Hpop; i++)
		fscanf(fp, "%d %d %u %u %u %u %f %d", &(HL+i)->X.x, &(HL+i)->X.y, &(HL+i)->EN, &(HL+i)->SN, 
														&(HL+i)->HB, &(HL+i)->DR, &(HL+i)->energy, &(HL+i)->age );

	for(int i=0; i<Cpop; i++)
		fscanf(fp, "%d %d %u %u %u %u %f %d", &(CL+i)->X.x, &(CL+i)->X.y, &(CL+i)->EN, &(CL+i)->SN, 
													&(CL+i)->CN, &(CL+i)->DR, &(CL+i)->energy, &(CL+i)->age );
	fclose( fp );
	
	for(int i=0; i<Ppop; i++)						// Link organisms into ecosystem grid
		(E +SX*(PL+i)->X.y + (PL+i)->X.x)->P = (PL+i);
		
	for(int i=0; i<Hpop; i++)
		(E +SX*(HL+i)->X.y + (HL+i)->X.x)->H = (HL+i);
		
	for(int i=0; i<Cpop; i++)
		(E +SX*(CL+i)->X.y + (CL+i)->X.x)->C = (CL+i);
	}

void storeSamples( void )			// Store 1000 samples of the genetics of P, H, C population
	{
	FILE *fp;
	fp = fopen( "sample.dat", "w" );
	for(int i=0; i<1000; i++)
		{
		int R = rand()%Ppop;			// Randomly select from the plant population
		fprintf(fp, "%u %u %u %f %d\n", (PL+R)->EN, (PL+R)->PL, (PL+R)->GM,
											  (PL+R)->energy, (PL+R)->age );
		}
	for(int i=0; i<1000; i++)
		{
		int R = rand()%Hpop;			// Randomly select from the herbivore population
		fprintf(fp, "%u %u %u %u %f %d\n", (HL+R)->EN, (HL+R)->SN, (HL+R)->HB, (HL+R)->DR,
											  		  (HL+R)->energy, (HL+R)->age );
		}		
	for(int i=0; i<1000; i++)
		{
		int R = rand()%Cpop;			// Randomly select from the carnivore population
		fprintf(fp, "%u %u %u %u %f %d\n", (CL+R)->EN, (CL+R)->SN, (CL+R)->CN, (CL+R)->DR,
													  (CL+R)->energy, (CL+R)->age );
		}
	fclose( fp );	
	}

	// Sample buffers
struct pnode *PS = (struct pnode *) malloc( 1000 * sizeof( struct pnode ));
struct hnode *HS = (struct hnode *) malloc( 1000 * sizeof( struct hnode ));
struct cnode *CS = (struct cnode *) malloc( 1000 * sizeof( struct cnode ));
	// Second sample buffers
struct pnode *PB = (struct pnode *) malloc( 1000 * sizeof( struct pnode ));
struct hnode *HB = (struct hnode *) malloc( 1000 * sizeof( struct hnode ));
struct cnode *CB = (struct cnode *) malloc( 1000 * sizeof( struct cnode ));

void loadSamples( void )
	{
	FILE *fp;
	fp = fopen( "sample.dat", "r" );
	
	for(int i=0; i<1000; i++)
		{
		fscanf(fp, "%u %u %u %f %d", &(PS+i)->EN, &(PS+i)->PL, &(PS+i)->GM, 
										  	  &(PS+i)->energy, &(PS+i)->age );
		if( fRan() < 0.27 )			// 7% marking rate for testing purposes only :)
			(PS+i)->GM = 3855;						// Maternal genetic marker
		else
			(PS+i)->GM = 0;
		}
			
	for(int i=0; i<1000; i++)
		fscanf(fp, "%u %u %u %u %f %d", &(HS+i)->EN, &(HS+i)->SN, &(HS+i)->HB, &(HS+i)->DR, 
												  &(HS+i)->energy, &(HS+i)->age );
				
	for(int i=0; i<1000; i++)
		fscanf(fp, "%u %u %u %u %f %d", &(CS+i)->EN, &(CS+i)->SN, &(CS+i)->CN, &(CS+i)->DR, 
												  &(CS+i)->energy, &(CS+i)->age );	
	
	fclose( fp );
	}

void loadSampTwo( void )							// Fill second sample buffers
	{
	FILE *fp;
	fp = fopen( "sample9.dat", "r" );			// Naming scheme ??   Hmm...
	
	for(int i=0; i<1000; i++)
		fscanf(fp, "%u %u %u %f %d", &(PB+i)->EN, &(PB+i)->PL, &(PB+i)->GM, 
										  	  &(PB+i)->energy, &(PB+i)->age );
			
	for(int i=0; i<1000; i++)
		fscanf(fp, "%u %u %u %u %f %d", &(HB+i)->EN, &(HB+i)->SN, &(HB+i)->HB, &(HB+i)->DR, 
												  &(HB+i)->energy, &(HB+i)->age );
				
	for(int i=0; i<1000; i++)
		fscanf(fp, "%u %u %u %u %f %d", &(CB+i)->EN, &(CB+i)->SN, &(CB+i)->CN, &(CB+i)->DR, 
												  &(CB+i)->energy, &(CB+i)->age );	
	
	fclose( fp );
	}
	
	// Create a new ecosystem using sampled genetic data.
void useSamples( void )
	{
	int X, Y;
		
	initEcosystem();									// Refresh the nutrient layer at initial levels

//	Ppop = 800000; Hpop = 200000; Cpop = 30000;	// Numbers from current run time
	Ppop = 500000; Hpop = 230000; Cpop = 50000;

	for(int i=0; i<Ppop/2; i++)					// Fill half the space from sample.dat
		{
		int R = rand()%1000;							// Choose a sampled plant at random
															// Sample files are 1000 units per organism
		do	{			// Unless there is a living plant in that location
			X = rand()%SZ;								// Generate new point
			Y = rand()%SZ; 
			} 	while( (E +SZ*Y +X)->P );			// When a plant is resident
		
		(PL+i)->X.x = X;
		(PL+i)->X.y = Y;
		(E +SZ*Y +X)->P = (PL+i);
		(PL+i)->EN = (PS+R)->EN;
		(PL+i)->PL = (PS+R)->PL;					// (PL+i)->GM = 3855;
		(PL+i)->GM = 3855; //(PS+R)->GM;			// Maternal genetic marker TODO	
															// From 'sample.dat'
		(PL+i)->energy = (PS+R)->energy;
		(PL+i)->age = (PS+R)->age;
		(PL+i)->alive = true;		
		}
	
	for(int i=Ppop/2; i<Ppop; i++)				// Use second samples
		{													
		int R = rand()%1000;							// Choose a sampled plant at random
		
		do	{			// Unless there is a living plant in that location
			X = rand()%SZ;								// Generate new point
			Y = rand()%SZ; 
			} 	while( (E +SZ*Y +X)->P );			// When a plant is resident
		
		(PL+i)->X.x = X;
		(PL+i)->X.y = Y;
		(E +SZ*Y +X)->P = (PL+i);
		(PL+i)->EN = (PB+R)->EN;
		(PL+i)->PL = (PB+R)->PL;					// (PL+i)->GM = 3855;
		(PL+i)->GM = 0; //(PB+R)->GM;				// Maternal genetic marker TODO		
															// From 'sample9.dat'
		(PL+i)->energy = (PB+R)->energy;
		(PL+i)->age = (PB+R)->age;
		(PL+i)->alive = true;		
		}
	
	for(int i=0; i<Hpop/2; i++)
		{
		int R = rand()%1000;
		
		do	{		// Unless there is a herbivore living in that location
			X = rand()%SZ;								// Generate new point
			Y = rand()%SZ; 
			} 	while( (E +SZ*Y +X)->H );
		
		(HL+i)->X.x = X;
		(HL+i)->X.y = Y;
		(E +SZ*Y +X)->H = (HL+i);
		(HL+i)->EN = (HS+R)->EN;
		(HL+i)->SN = (HS+R)->SN;
		(HL+i)->HB = (HS+R)->HB;
		(HL+i)->DR = (HS+R)->DR;
		(HL+i)->energy = (HS+R)->energy;
		(HL+i)->age = (HS+R)->age;
		(HL+i)->alive = true;		
		}
	
	for(int i=Hpop/2; i<Hpop; i++)
		{
		int R = rand()%1000;
		
		do	{		// Unless there is a herbivore living in that location
			X = rand()%SZ;								// Generate new point
			Y = rand()%SZ; 
			} 	while( (E +SZ*Y +X)->H );
		
		(HL+i)->X.x = X;
		(HL+i)->X.y = Y;
		(E +SZ*Y +X)->H = (HL+i);
		(HL+i)->EN = (HB+R)->EN;
		(HL+i)->SN = (HB+R)->SN;
		(HL+i)->HB = (HB+R)->HB;
		(HL+i)->DR = (HB+R)->DR;
		(HL+i)->energy = (HB+R)->energy;
		(HL+i)->age = (HB+R)->age;
		(HL+i)->alive = true;		
		}
		
	for(int i=0; i<Cpop/2; i++)
		{
		int R = rand()%1000;
		
		do	{		// Unless there is a carnivore living in that location
			X = rand()%SZ;								// Generate new point
			Y = rand()%SZ; 
			} 	while( (E +SZ*Y +X)->C );
		
		(CL+i)->X.x = X;
		(CL+i)->X.y = Y;
		(E +SZ*Y +X)->C = (CL+i);
		(CL+i)->EN = (CS+R)->EN;
		(CL+i)->SN = (CS+R)->SN;
		(CL+i)->CN = (CS+R)->CN;
		(CL+i)->DR = (CS+R)->DR;
		(CL+i)->energy = (CS+R)->energy;
		(CL+i)->age = (CS+R)->age;
		(CL+i)->alive = true;
		}
		
	for(int i=Cpop/2; i<Cpop; i++)
		{
		int R = rand()%1000;
		
		do	{		// Unless there is a carnivore living in that location
			X = rand()%SZ;								// Generate new point
			Y = rand()%SZ; 
			} 	while( (E +SZ*Y +X)->C );
		
		(CL+i)->X.x = X;
		(CL+i)->X.y = Y;
		(E +SZ*Y +X)->C = (CL+i);
		(CL+i)->EN = (CB+R)->EN;
		(CL+i)->SN = (CB+R)->SN;
		(CL+i)->CN = (CB+R)->CN;
		(CL+i)->DR = (CB+R)->DR;
		(CL+i)->energy = (CB+R)->energy;
		(CL+i)->age = (CB+R)->age;
		(CL+i)->alive = true;
		}		
	}

void showNew( void )
	{
	for(int i=Ppop-10; i<Ppop; i++)	
		printf("%u %u %u\n", (PL+i)->EN, (PL+i)->PL, (PL+i)->GM );
	}

void showSamples( void )
	{
	chromosome EN, PL, GM, CN, HB, SN, DR;
	float energy;
	int age;
	FILE *fp;
	fp = fopen( "sample.dat", "r");
	
	for(int i=0; i<1000; i++)
		{
		fscanf(fp, "%u %u %u %f %d", &EN, &PL, &GM, &energy, &age );
		geneLog( EN );	printf("%s\n", G );
		geneLog( PL );	printf("%s\n", G );
		geneLog( GM );	printf("%s\n", G );
		}
	printf("\n");
		
	for(int i=0; i<1000; i++)
		{
		fscanf(fp, "%u %u %u %u %f %d", &EN, &SN, &HB, &DR, &energy, &age );
		geneLog( EN );	printf("%s\n", G );
		geneLog( SN );	printf("%s\n", G );
		geneLog( HB );	printf("%s\n", G );
		geneLog( DR );	printf("%s\n", G );
		}
	printf("\n");
		
	for(int i=0; i<1000; i++)
		{
		fscanf(fp, "%u %u %u %u %f %d", &EN, &SN, &CN, &DR, &energy, &age );
		geneLog( EN );	printf("%s\n", G );
		geneLog( SN );	printf("%s\n", G );
		geneLog( CN );	printf("%s\n", G );
		geneLog( DR );	printf("%s\n", G );
		}
		
	fclose( fp );
	}

void stamp( void )
	{
	time_t t;
   struct tm *nt;

   time( &t );
   nt = localtime( &t );
   
   if( nt->tm_hour < 10 )
   	printf("0");
   printf("%d:", nt->tm_hour );
   
   if( nt->tm_min < 10 )
   	printf("0");
   printf("%d:", nt->tm_min );
   
   if( nt->tm_sec < 10 )
   	printf("0");
   printf("%d\n", nt->tm_sec );
	}

void fStamp( FILE *fp )
	{
	time_t t;
	struct tm *nt;
	
	time( &t );
	nt = localtime( &t );
	
	if( nt->tm_hour < 10 )
   	fprintf( fp, "0");
   fprintf( fp, "%d:", nt->tm_hour );
   
   if( nt->tm_min < 10 )
   	fprintf( fp, "0");
   fprintf( fp, "%d:", nt->tm_min );
   
   if( nt->tm_sec < 10 )
   	fprintf(fp, "0");
   fprintf( fp, "%d\n", nt->tm_sec );
	}

/////////////////////////////////////////////////////////////////////////////////////////////////
///////////////  Rarely Used Functions  /////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////

void mergeSamples( void )
	{
	int X, Y;
	// loadData() reads eco.dat and creates the ecosystem
	// Open the sample.dat file
	// Add those organisms to the pre-existing ecosytem
	// Mark the plant chromosome GM of the sample data.
	
	printf("%d %d %d\n", Ppop, Hpop, Cpop );
	
	for(int i=0; i<990000; i++)						// Add plants from sample.dat
		{
		int R = rand()%1000;							// Choose a sampled plant at random
//		printf("R = %d\n", R );
		int j = Ppop + i;								// Offset to end of plant list
		
		do	{			// Unless there is a living plant in that location
			X = rand()%SZ;								// Generate new point
			Y = rand()%SZ; 
			} 	while( (E +SZ*Y +X)->P );			// When a plant is resident
		
		(PL+j)->X.x = X;		// Ppop+i
		(PL+j)->X.y = Y;
		(E +SZ*Y +X)->P = (PL+j);
		(PL+j)->EN = (PS+R)->EN;
		(PL+j)->PL = (PS+R)->PL;					// (PL+i)->GM = 3855;
		(PL+j)->GM = 3855;							// Maternal genetic marker 		
		
		(PL+j)->energy = (PS+R)->energy;
		(PL+j)->age = (PS+R)->age;
		(PL+j)->alive = true;
		}
	Ppop += 990000;										// Update plant population count
	
	for(int i=0; i<1000; i++)
		{
		int R = rand()%1000;
		int j = Hpop + i;								// Offset to end of herbivore list
		
		do	{		// Unless there is a herbivore living in that location
			X = rand()%SZ;								// Generate new point
			Y = rand()%SZ; 
			} 	while( (E +SZ*Y +X)->H );
		
		(HL+j)->X.x = X;
		(HL+j)->X.y = Y;
		(E +SZ*Y +X)->H = (HL+j);
		(HL+j)->EN = (HS+R)->EN;
		(HL+j)->SN = (HS+R)->SN;
		(HL+j)->HB = (HS+R)->HB;
		(HL+j)->DR = (HS+R)->DR;
		(HL+j)->energy = (HS+R)->energy;
		(HL+j)->age = (HS+R)->age;
		(HL+j)->alive = true;		
		}
	Hpop += 1000;
		
	for(int i=0; i<1000; i++)
		{
		int R = rand()%1000;
		int j = Cpop + i;
		
		do	{		// Unless there is a carnivore living in that location
			X = rand()%SZ;								// Generate new point
			Y = rand()%SZ; 
			} 	while( (E +SZ*Y +X)->C );
		
		(CL+j)->X.x = X;
		(CL+j)->X.y = Y;
		(E +SZ*Y +X)->C = (CL+j);
		(CL+j)->EN = (CS+R)->EN;
		(CL+j)->SN = (CS+R)->SN;
		(CL+j)->CN = (CS+R)->CN;
		(CL+j)->DR = (CS+R)->DR;
		(CL+j)->energy = (CS+R)->energy;
		(CL+j)->age = (CS+R)->age;
		(CL+j)->alive = true;
		}
	Cpop += 1000;
	printf("%d %d %d\n", Ppop, Hpop, Cpop );
	}

void display( chromosome seq )
	{
	int j=0;
	chromosome mask = 1 << 31;						// Shift mask all the way left
	chromosome temp = seq;
	for(int i=0; i<32; i++)
		{
		if( seq & mask )								// Compare mask bit to present bit
			printf("1");
		else
			printf("0");
		if( j++ == 3 )									// Separate into bit quartets
			{
			printf(" ");
			j=0;
			}
		seq <<= 1;										// Shift sequence to the left
		}
	printf("\n");
	}

void geneLog( chromosome seq )					// Log sample genes to disk
	{
	unsigned int mask = 1 << 31;					// Shift mask all the way left

	for(int i=0; i<32; i++)
		{
		if( seq & mask )								// Compare mask bit to present bit
			G[i] = '1';
		else
			G[i] = '0';

		seq <<= 1;										// Shift sequence to the left
		}
	}

void showGenes( void )								// RAND_MAX = 2147483647
	{	
	FILE *fp;
	fp = fopen("bin.dat","w");
	
		// Choose 10 random plants with rand()%Ppop
	for(int i=0; i<10; i++)							// Instant karma's gonna get you
		{													// Gonna knock you right on the head
		int R = rand()%Ppop;							// You better get yourself together
															// Pretty soon you're gonna be dead
		geneLog( (PL+R)->EN );
		fprintf( fp,"%s\n", G );
		
		geneLog( (PL+R)->PL );
		fprintf( fp,"%s\n", G );
		}
		
	for(int i=0; i<10; i++)
		{
		int R = rand()%Hpop;
		
		geneLog( (HL+R)->EN );
		fprintf( fp,"%s\n", G );
		geneLog( (HL+R)->SN );
		fprintf( fp,"%s\n", G );
		geneLog( (HL+R)->HB );
		fprintf( fp,"%s\n", G );
		geneLog( (HL+R)->DR );
		fprintf( fp,"%s\n", G );
		}

	for(int i=0; i<10; i++)				
		{	
		int R = rand()%Cpop;	
		
		geneLog( (CL+R)->EN );
		fprintf( fp,"%s\n", G );
		geneLog( (CL+R)->SN );
		fprintf( fp,"%s\n", G );
		geneLog( (CL+R)->CN );
		fprintf( fp,"%s\n", G );
		geneLog( (CL+R)->DR );
		fprintf( fp,"%s\n", G );
		}	
	fclose( fp );
	}

//	showGenes();										// Store genetic samples in 'bin.dat'
//	char c;												// Dummy character buffer										
//	scanf("%c", &c );									// Gatekeeper :)

void checkEnergy( void )	
	{
	printf("%d %d %d\n", Ppop, Hpop, Cpop );	
//	for(int i=0; i<Ppop; i++)
//		printf("%d  %f\n", i, (PL+i)->energy );
	for(int i=0; i<Hpop; i++)
		if( (HL+i)->energy > 0.0 ) printf("%d  %f\n", i, (HL+i)->energy );
//	for(int i=0; i<Cpop; i++)
//		printf("%d  %f\n", i, (CL+i)->energy );
	} 
	
void checkOffspring( void )
	{
	printf("%d %d %d\n", newPpop, newHpop, newCpop );	
	for(int i=0; i<newPpop; i++)
		printf("P%d  %d %d\n", i, (NP+i)->X.x, (NP+i)->X.y );
	for(int i=0; i<newHpop; i++)
		printf("H%d  %d %d\n", i, (NH+i)->X.x, (NH+i)->X.y );
	for(int i=0; i<newCpop; i++)
		printf("C%d  %d %d\n", i, (NC+i)->X.x, (NC+i)->X.y );
	}

void checkOthers( void )
	{
	for(int i=0; i<Ppop; i++)
		if( i % 10 == 0 )								// Decimate list
			printf("p %6.2f  ", (OPL+i)->energy );
	printf("\n");
	
	for(int i=0; i<Hpop; i++)
		if( i % 10 == 0 )
			printf("h %6.2f  ", (OHL+i)->energy );
	printf("\n");
	
	for(int i=0; i<Cpop; i++)
		if( i % 10 == 0 )
			printf("c %6.2f  ", (OCL+i)->energy );
	printf("\n");	
	}

void checkAlternate( void )						// Sample new organism lists
	{
	for(int i=0; i<Ppop; i++)
		if( (i % 50) == 0 )
			printf("%4d  P %4d %4d %9.2f\n", i, (OPL+i)->X.x, (OPL+i)->X.y, (OPL+i)->energy );
			
	for(int i=0; i<Hpop; i++)
		if( (i % 50) == 0 )
			printf("%4d  H %4d %4d %9.2f\n", i, (OHL+i)->X.x, (OHL+i)->X.y, (OHL+i)->energy );
			
	for(int i=0; i<Cpop; i++)
		if( (i % 50) == 0 )
			printf("%4d  C %4d %4d %9.2f\n", i, (OCL+i)->X.x, (OCL+i)->X.y, (OCL+i)->energy );
	}

void checkDir( void )
	{
	for(int i=0; i<9; i++)							// Check dir[]
		printf("%d  %d %d\n", i, dir[i].x, dir[i].y );
	}
	
void displayNew( void )								// Sample new organism lists
	{
	for(int i=0; i<newPpop; i++)
		if( (i % 50) == 0 )
			printf("%4d %9.2f\n", i, (NP+i)->energy );
	for(int i=0; i<newHpop; i++)
		if( (i % 50) == 0 )
			printf("%4d %9.2f\n", i, (NH+i)->energy );
	for(int i=0; i<newCpop; i++)
		if( (i % 50) == 0 )
			printf("%4d %9.2f\n", i, (NC+i)->energy );
	}

void check( void )									// Test SZ is 500 => 250,000 
	{
	for(int i=0; i<SZ*SZ; i++)						// Display organisms in Ecosystem
		{
		struct enode *Q = (E+i);					// Organism energy gene and levels
		if( Q->P && Q->P->alive )
			printf("E %6d  P %4d %4d %7.2f %7.2f\n", i, Q->P->X.x, Q->P->X.y, 
														(float) (Q->P->EN & 7) + 6.1, Q->P->energy );		
		if( Q->H && Q->H->alive )
			printf("E %6d  H %4d %4d %7.2f %7.2f\n", i, Q->H->X.x, Q->H->X.y,
			 											(float) (Q->H->EN & 63) + 20.2, Q->H->energy );		
		if( Q->C && Q->C->alive )
			printf("E %6d  C %4d %4d %7.2f %7.2f\n", i, Q->C->X.x, Q->C->X.y,
			 											(float) (Q->C->EN & 127) + 15.3, Q->C->energy );
		}
	}

bool chck( void )										// Test SZ is 500 => 250,000 
	{
	bool alive = false;
	for(int i=0; i<SZ*SZ; i++)						// Display organisms in Ecosystem
		{
		struct enode *Q = (E+i);
		if( Q->P && Q->P->alive )	alive = true;	
		if( Q->H && Q->H->alive )	alive = true;
		if( Q->C && Q->C->alive )	alive = true;
		}
	return !alive;
	}

void printPlants( void )
	{
	for(int j=0; j<SZ; j++)
		{
		for(int i=0; i<SZ; i++)
			{
			if( (E + j*SZ + i)->P )
				printf("P");
			else
				printf(".");
			}
		printf("\n");
		}	
	}

void printHerbivores( void )
	{
	for(int j=0; j<SZ; j++)
		{
		for(int i=0; i<SZ; i++)
			{
			if( (E +j*SZ +i)->H )
				printf("H");
			else
				printf(".");
			}
		printf("\n");
		}	
	printf("\n");
	}
	
void printCarnivores( void )
	{
	for(int j=0; j<SZ; j++)
		{
		for(int i=0; i<SZ; i++)
			{
			if( (E +j*SZ +i)->C )
				printf("C");
			else
				printf(".");
			}
		printf("\n");
		}	
	printf("\n");
	}

/////////////////////////////////////////////////////////////////////////////////////////////////
///////////////  Notes  /////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////

/*  From move()

//int hist[9];											// Histogram array

	static int first = 0;
	if( first == 0 )									// Display probabilities on first pass only
		{
		for(int i=1; i<9; i++)
			printf("%7.4f  ", dirProb[i] );
		printf("\n");
		first++;
		}

//	printf("%7.1f  ", tt );							// Probability total
//	printf("%7.4f  ", rn );							// Probability choice
	printf("%2d ", dr );								// Direction choice
	
	switch( dr )
		{
		case 1:  printf("East      ");	break;
		case 2:  printf("West      ");	break;
		case 3:  printf("North     ");	break;
		case 4:  printf("South     ");	break;
		case 5:  printf("Northeast ");	break;
		case 6:  printf("Southeast ");	break;
		case 7:  printf("Northwest ");	break;
		case 8:  printf("Southwest ");	break;
		}
	printf("\n");
*/				
//	hist[ dr ]++;										// Create a histogram of directions

// When you add graphics to this application
// Pick a good background color
// and three organism colors
// Then, add a brighter second color for each organism
// This signifies the newest plant, herbivore, or carnivore
// respectively.  Fade as age increases to standard organism color
// Make sure to layer the trophic level displays
// One each for nutrient, plant, herbivore, and carnivore levels.
// Then you can use 3D graphics, with rotations/translations or course :)
// to pick and display the layers you want to view.
// You could even pick different traits for each layer.  
// Energy or age come to mind. 

#ifdef NOTES
Missing Mustelidae:

Erascible Ermine
Fighting Fishers
Frenzied Ferrets
Mighty Mink
Marauding Martens
Oscillating Otters
Stout State Stoats
Warrior Weasels
#endif

// Just power down system first, install drive, and Linux should recognize it upon boot. 
// Then you just need to format it, mount it, and (optionally) 
// add the mount to fstab to have it persist across reboots.

// sudo dd if=/dev/sdc of="/media/kevin/H Back Two/BACK/Hermes.img"

//	RAND_MAX = 2147483647;

//		printf(".");		Progress markers
//		printf("*");
//		printf("|");
//		printf("^");
//		printf("+");
//		printf("\n");

#ifdef NOTES
void eat( void )
	{							
	struct pnode *P;									// Point at plants
	struct hnode *H;									// Point at herbivores
	struct vt A, B;									// Working vectors
	int X, Y;											// Location indices
	int j, ring;										// Direction, ring#
	float satiety, absorb;							// Genetic expression buffers	
	
		// Herbivores eat at E(x,y) and the five rings around it
	for(int i=0; i<Hpop; i++)
		{
		j = 0;											// Don't skip center cell
		ring = 1;
		A = (HL+i)->X;									// Current location		
		while( true )									// Forever loop needs a break out
			{												// Find an open spot within one cell			
			B = scale( ring, dir[j] );				// Motion vector
			X = A.x + B.x;								// Prospective new location 
			Y = A.y + B.y;

			if( X < 0 ) 			X += SZ;			// Wrap x,y indices
			else if( X > SZ-1 ) 	X -= SZ;			// into a toroidal space
			if( Y < 0 ) 			Y += SZ;
			else if( Y > SZ-1 ) 	Y -= SZ;

			satiety = (float) (((HL+i)->EN & (15 << 10)) >> 10) + 5.0; // EN::G3	5.0 - 20.0			
			if( (E +Y*SZ +X)->P && ((HL+i)->energy < satiety ))
				{
				P = (E +Y*SZ +X)->P;  				// Point to plant at E(x,y)	
				absorb = (float) ((((HL+i)->EN & (7 << 7)) >> 7) + 1.0) / 13.0; // EN::G4  	7.7 - 61.5%	
				printf("h");
				(HL+i)->energy += P->energy * absorb;	// Add energy to herbivore
				P->energy -= P->energy * absorb;			// Remove energy from plant
				(E +Y*SZ +X)->nutrient += P->energy;	// Remainder to nutrient layer
				(E +Y*SZ +X)->P = NULL;				// Plant no longer exists here
				P->alive = false;						// Plant is deceased		
				}
						
			j++;											// Walk around direction vector				
			if( j > 8 ) 								// Limit search count
				{											// 8 random choices / ring
				j = 1;									// Search next
				ring++;									// outer ring
				}
			if( ring > 5 ) break;					// PFA step limiter  NEW			
			}			
		}
	
		// Carnivores eat in rings around E(x,y) too
	for(int i=0; i<Cpop; i++)
		{
		j = 0;											// Don't skip center cell
		ring = 1;
		A = (CL+i)->X;
		while( true )									// Eat until break at ring 12
			{
			B = scale( ring, dir[j] );				// Motion vector
			X = A.x + B.x;
			Y = A.y + B.y;
			
			if( X<0 ) 				X += SZ;			// Toroidal space
			else if( X>SZ-1 ) 	X -= SZ;
			if( Y<0 ) 				Y += SZ;
			else if( Y>SZ-1 ) 	Y -= SZ;
			
			satiety = (float) (((CL+i)->EN & (15 << 10)) >> 10) + 8.0; // 8 - 23.0
			if( (E +Y*SZ +X)->H && ((CL+i)->energy < satiety ))
				{				
				H = (E +Y*SZ +X)->H;					// Point to herbivore at E(x,y)	
				absorb = (float) ((((CL+i)->EN & (7 << 7)) >> 7) + 1.0) / 18.0;	
				printf("c");	
				(CL+i)->energy += H->energy * absorb;	// Add energy to herbivore
				H->energy -= H->energy * absorb;			// Remove energy from plant				
				(E +Y*SZ +X)->nutrient += H->energy;	// Remainder to nutrient layer					
				(E +Y*SZ +X)->H = NULL;				// Herbivore no longer exists here					
				H->alive = false;						// Herbivore is deceased								
				}
		
			j++;											// Walk around direction vector				
			if( j > 8 ) 								// Limit search count
				{											// 8 random choices / ring
				j = 1;									// Search next
				ring++;									// outer ring
				}
			if( ring > 9 ) break;					// PFA step limiter
			}
		}
	}
#endif
#ifdef NOTES		
Diffuse deposited nutrients over time:
	Nutrients are supplied through replenishment, feeding, or reproduction.
	Need nutrient diffusion time step relationship
	j = - D * dc/dx
	D is the diffusivity constant
	dc/dx is rate of change of concentration in the x direction (negative away )
	And j is the amount flowing through a space
	
	How about a little variation in nutrient replenishment?
	Richer areas along watersheds, or in patches
	Add along one or two edges with rhythmic pulses
	Flow from one or two edges toward and out of the opposite side or sides
	Flow could be toroidal or not, depending on your whims :)
#endif

//		char c;											// Dummy character buffer
//		scanf("%c", &c );								// Gatekeeper :)
//		if( kk%10 == 0 ) scanf("%c", &c );		// Gate keeper, Keymaster :)

//		if( kk%100 == 0 )								// Store gene samples every 100 cycles
//			showGenes();								// In 'bin.dat' 10 plant, 10 herb, and 10 carn

#ifdef NOTES
This would be called with: (CL+i)->Hsense[]  
It would scan the sense array for the highest number
and return its index.  Use that as the index to dir[ k ]
and use scale and ring and rand()%speed to move to that spot.
That would get the carnivore moving into the most dense 
area of the nearby herd.
Remember to avoid leaping onto a fellow carnivore.  It is rude.

void carnScan( void )								// Sense herbivores in local area
	{														// 8*8 = 64 possible cells
	int maxRing = 8;									// PFA number	potential gene  111 + 1								
	struct vt A, B, C;								// Working vectors
	
	for(int i=0; i<Cpop; i++)
		{
		A = (CL+i)->X;									// You are here.
										
		for(int k=1; k<9; k++)						// Clear Hsense[]
			(CL+i)->Hsense[ k ] = 0;
			
		for(int ring=1; ring<maxRing; ring++)
			for(int k=1; k<9; k++)
				{
				B = scale( ring, dir[ k ] );		// Walk out the direction vector
				C = bound( sum( A, B ) );			// Search location				
				if( (E +C.y*SZ +C.x)->H )  (CL+i)->Hsense[ k ]++;	// I found a herbivore :)
				}
		}
	}

// hiLo.c
// K. J. Rock
// October 28, 2024

// g++ -o HL hiLo.c

#include <time.h>
#include <stdio.h>
#include <stdlib.h>

// I think (HL+i)->Psense[] and (CL+i)->Hsense[] are sufficient.
// Then it would be feeding and pouncing with no running.
// That way there is no need a higher metabolic rate
// Then use maxPlant() and maxHerb()
// How do I find the second highest density paths?

int maxHerb( int CS[] )
	{
	int Hmax;
	int max = -100;		// Guard
	for(int i=1; i<9; i++)
		if( CS[ i ] > max )
			{
			Hmax = i;
			max = CS[ Hmax ]; 
			}
	return Hmax;
	}

int minHerb( int HS[] )								// Should really be
	{														// minCarn( ) to avoid them
	int Hmin;
	int min = 100;			// Guard
	for(int i=1; i<9; i++)
		if( HS[ i ] < min )
			{
			Hmin = i;
			min = HS[ Hmin ];
			}
	return Hmin;
	}

int main( void )
	{	
	int Hmax, Hmin;	
	int HS[9];											// Herbivore Sensor array
	
	srand( (unsigned) time( NULL ) );			// Seed random number generator
	
	for(int k=0; k<21; k++)							// Loop for effect
		{
		for(int i=0; i<9; i++)
			HS[ i ] = rand()%15;						// Fill array with random values
		
		Hmin = minHerb( HS );						// Find the least dense direction
		printf("HS[%1d] = %1d\t", Hmin, HS[ Hmin ] );
		
		Hmax = maxHerb( HS );						// Find the most dense direction																					
		printf("HS[%1d] = %2d\n", Hmax, HS[ Hmax ] );	// Display maximum 
		}
	}
#endif
#ifdef NOTES
			Five senses { sight, hearing, smell, taste, touch }
Sensor genes:      ->SN  in herbivores and carnivores
  	8 bits sight				camouflage effectiveness would be tested here
	8 bits hearing				hear predator or prey
	8 bits smell     32-24 = 8		smell poison or prey
	4 bits taste				taste poison, touch spines
	4 bits touch		2^4 = 16  pursue vibration of ground?? or hide if sensed 

see =   (float)  ((CL+i)->SN & 255) / 255.0;						SN::G5  sight range
hear =  (float) (((CL+i)->SN & (255 << 8)) >> 8) / 255.0;	SN::G4  hearing ability
smell = (float) (((HL+i)->SN & (255 << 16)) >> 16) / 255.0;	SN::G3  scent capability
taste	= (float) (((HL+i)->SN & (15 << 24)) >> 24) / 15.0;	SN::G2  taste sense
touch = (float) (((HL+i)->SN & (15 << 28)) >> 28) / 15.0;	SN::G1  touch sensitivity

// Hard wire plant senescence at 5 cycles
// Then set its maturity at 3 or 4 cycles.
// Thus each could create two or three offspring.
// They can get eaten at any age from 1 to 5 cycles.

pollen = (float) ((PL+i)->PL & 63) / 63;                    PL::G5  pollination probability
#endif

#ifdef NOTES		
void moveHerbivores( void )
	{
	struct vt A, B;									// Working vectors
	int X, Y;
	int min, ring, ch, speed;

	for(int i=0; i<Hpop; i++)						// Scan herbivore list
		{
		ch = 0;											// Limit choice count
		min = 1;											// 
															// { 1 <= speed <= 8 }
		speed = (int) (((HL+i)->HB & (7 << 28)) >> 28 ) + min;  // 7 + 29 = 36
		
		// Need to express (HL+i)->Psense[] & (HL+i)->Csense[]
		// to determine next move.  Hmm.. ??
		// We could scan our sense arrays in order
		// If Xsense[ ? ] is greater than zero 
		// Hold direction steady and step ring until maxRing limit
		// Then step direction to the next non-zero choice
		// Or, conversely, we can avoid carnivores by moving
		// in the least populated directions.
		
// Sensor scan algorithm.
// You will need  a version for herbivores
// and another for carnivores since they don't seek the same goals.

If a carnivore senses a high density of prey
move will respond by moving in that direction

If a herbivore senses a carnivore nearby it enables
the high metabolic rate gene and increases its speed.
Move() will aim in the lowest carnivore density direction
at the new rate of speed.  Move() also extracts the added
energy required by running.  Now, how long do I run away?
Ah, there is a gene for that :)
		
			// Scan from inner ring to outer ring defined by sensor strength
		for(int ring=1; ring<maxRing; ring++)			
			for(int k=1; k<9; k++)					// Scan directions
				{				
				B = scale( ring, dir[ k ] );  	// Walk around ring 1..8
				C = bound( sum( A, B ) );			// Search location

				if( (E +C.y*SZ +C.x)->P)  (HL+i)->Psense[ k ]++;	// I found a plant
				if( (E +C.y*SZ +C.x)->C)  (HL+i)->Csense[ k ]++;	// I found a carnivore
				}
		
Remember a herbivore eats at dir[ 0 ] too.
Start at ring = 0 for this to be true

Walk around (HL+i)->Psense[ k ] > 0 
Use that direction choice 
B = scale( ring, dir[ k ] );
Step ring++ to maxRing  or value of (HL+i)->Psense[] use a counter?

Repeat with (HL+i)->Csense[ k ] > 0 to override prospective new location.
		
	
		A = (HL+i)->X;									// Current location		
		do	{												// Find an open spot within range	
			ring = min + rand()%speed;				// {min <= ring < speed}		
			B = scale( ring, move( (HL+i)->DR ) );	// Express direction genes
			C = bound( sum( A, B ) );			// Search location
			
			if( ch > 8 ) break; else ch++;		// Limit search count to 8 choices
			} 	while( (E +C.y*SZ +C.x)->H );			// While location is occupied
		
		(E +A.y*SZ +A.x)->H = NULL;				// Empty old location
		(HL+i)->X.x = X;								// Remember where I live
		(HL+i)->X.y = Y;								// Index into Ecosystem coordinates
		(E +C.y*SZ +C.x)->H = (HL+i);					// Fill new location			
		}
	}

In moveCarnivore()  :

Add pursuit speed here.  We could call the metabolism() function here
Then control resting, active, or chase metabolic rates.  Yes, Three of them.
How do I used the speed gene to move along a target rich path?
I think find the highest count and move in that direction will work
Remember to use a check for occupied cells.  Carnivores avoid carnivores

OK, a plan: determine highest ->Hsense[ k ] 
Then move along dir[ k ] using scale and ring.
That should work.

It would seek the highest local concentration of prey.
I think that would be a good time to switch from resting
to active metabolism.
We can try pursuit metabolism after this works.	
#endif
#ifdef NOTES
Add GM genetic marker chromosome 
GM = 111100001111 = 3855 or 10101010 = 170

Implement gene marker which mother will pass on.  
This chromosome will not run through crossover()
What should be in it?
Take samples every XX cycles
See how the marker is spread around those samples.
GM = 1111 0000 1111 0000 1111 0000 1111 0000 ;   for a distinctive pattern
11110000111100001111000011110000
NO, use a small number instead.
Such as GM = 111100001111 = 3855 or 10101010 = 170
Sample data shows chromosomes are all 10 digit numbers
3855 or 170 would stick out nicely.
#endif

// https://archive.is/s0Szb  MEMS GHz filters
// https://www.3blue1brown.com/lessons/holograms
// https://www.nobelprize.org/uploads/2018/06/gabor-lecture.pdf
// https://www.scientificamerican.com/article/a-cubic-millimeter-of-a-human-brain-has-been-mapped-in-spectacular-detail/
// https://bpb-us-e1.wpmucdn.com/websites.harvard.edu/dist/f/49/files/2022/08/H01_selectedneurons_mip4rothalfcrop2-2400x6060-1.png
// https://www.star.nesdis.noaa.gov/goes/conus_band.php?sat=G18&band=GEOCOLOR&length=48&dim=1
	
#ifdef NOTES
COLORREF SetPixel( bkhdc, X, Y, COLORREF color );

for(int j=0; j<SZ; j++)
	for(int i=0; i<SZ; i++)
		(E+j*SZ+i)->nutrient
		
		convert to color on temperature palette.
		Black body or Blue to Red to Yellow to White?
		SetPixel( bkhdc, i, j, color );

This will maintain a colored landscape representing the nutrient levels
We could modify the palette scaling by recording the 
energy minimum and maximum levels.  

void drawFrame( void )		// Draw frame of simulation
	{
	double ux, uy, uNorm, rho;

	for (int iY=2; iY<=sim.ly; ++iY) 
		for (int iX=1; iX<=sim.lx; ++iX) 
			{
			computeMacros( sim.lattice[iX][iY].fPop, &rho, &ux, &uy );
			uNorm = sqrt( ux*ux+uy*uy );
			
			// If this cell is not an obstruction of some sort.
			if( sim.lattice[iX][iY].ob == false )
				{
				if( vmax < uNorm ) vmax = uNorm;	// find velocity bounds
				if( vmin > uNorm ) vmin = uNorm;

				if( pmax < rho ) pmax = rho;		// find density bounds
				if( pmin > rho ) pmin = rho;
				}			
			}
	
	for (int iY=2; iY<=sim.ly; ++iY) 
		for (int iX=1; iX<=sim.lx; ++iX) 
			{			
			int x = 400 + iX;						// offset channel
			int y = 300 + iY;						// on display area
			
			int range = 650;						// hiColor
						
			int t;									// Color index
			t = colorTransform( range, 0, vmax, vmin, uNorm );	// velocity
		
			SetPixel( bkhdc, x, y, RGB( palette[t].R, palette[t].G, palette[t].B ) );
			}
	
//	drawTemplate();	
	}
	
// returns index into palette[i]
int colorTransform(int hiColor, int loColor, double hiProblem, double loProblem, double xw)
	{
	return((hiColor-loColor) / (hiProblem-loProblem) * (xw-loProblem) + loColor);
	}
	
void initPalette( void )
	{
	int q = 0;
	
	for(int i=0; i<256; i++)	// from blue to red
		{
		palette[q].R = i;
		palette[q].G = 0;
		palette[q].B = 255-i;
		q++;
		}
			
	for(int i=0; i<256; i++)	// from red to yellow
		{
		palette[q].R = 255;
		palette[q].G = i;
		palette[q].B = 0;
		q++;
		}
			
	for(int i=0; i<256; i++) 	// from yellow to white
		{
		palette[q].R = 255;
		palette[q].G =	255;
		palette[q].B = i;
		q++;
		}		// for 768 colors from blue to white
	}

	// Create a palette of 768 colors from black to white 
void initPalette( void )		// Black body radiation palette
	{
	int k = 0;											// Index 768 colored pens
	
		// create palette of pens
//	for (int i=0; i<256; i++)						// Blue to red transition
//		pn[k++] = CreatePen( PS_SOLID, 1, RGB( i, 0, 255-i) );
	
	for (int i=0; i<256; i++)						// Black to red transition
		pn[k++] = CreatePen( PS_SOLID, 1, RGB( i, 0, 0) );
		
	for (int i=0; i<256; i++)						// Red to yellow transition	
		pn[k++] = CreatePen( PS_SOLID, 1, RGB( 255, i, 0) );
	
	for (int i=0; i<256; i++)						// Yellow to white transition
		pn[k++] = CreatePen( PS_SOLID, 1, RGB( 255, 255, i) );
		
	int q = 0;											// Index brushes too.

		// create a palette of brushes
//	for (int i=0; i<256; i++)						// Blue to red
//		br[q++] = CreateSolidBrush( RGB( i, 0, 255-i) );
	
	for (int i=0; i<256; i++)						// Black to red transition
		br[q++] = CreateSolidBrush( RGB( i, 0, 0) );
		
	for (int i=0; i<256; i++)						// Red to yellow
		br[q++] = CreateSolidBrush( RGB( 255, i, 0) );
		
	for (int i=0; i<256; i++)						// Yellow to white
		br[q++] = CreateSolidBrush( RGB( 255, 255, i) );
	}		// k = q = 768 colors in B R Y W palette
	
#endif
#ifdef NOTES

When you add graphics to this application:
Pick a good background color and three organism colors.
Then, add a brighter second color for each organism.
This signifies the newest plant, herbivore, or carnivore respectively.  
Fade them as age increases, to standard organism color.
Make sure to layer the trophic level displays: 
One each - for nutrient, plant, herbivore, and carnivore levels.

Pick and display the layers you want to view.
You could even pick different traits for each layer.  
Energy or age come to mind.  
Track the energy stored at each point on every layer as a color map.

Adjust the color scale for each cycle

Build the black body radiation palette: Black to Red to Yellow to White

in replenish()
	{	// Determine energy extrema
	if( (E+i)->nutrient < emin ) emin = (E+i)->nutrient;
	if( (E+i)->nutrient > emax ) emax = (E+i)->nutrient;
	}

	// Fill nutrient layer bitmap
for(int i=0; i<SZ*SZ; i++)
	{
	t = colorTransform( range, 0, emax, emin, (E+i)->nutrient );	// Energy
	SetPixel( bkhdc, x, y, RGB( palette[t].R, palette[t].G, palette[t].B ) );
	}
	

hdc = GetDC( hwnd);		
bkhdc = CreateCompatibleDC( hdc);
bkBitmap = CreateCompatibleBitmap( hdc, WORLDX, WORLDY);
	Need a bitmap for nutrient layer
	and another bitmap for PHC layer

SelectObject( bkhdc, bkBitmap);     // Add a bitmap to the background hdc	
GetClientRect( hwnd, &rct);	
FillRect( bkhdc, &rct, (HBRUSH) (GetStockObject(GRAY_BRUSH)));


int color;											// Color index for pn[ 768 ]
	
	color = (int) 550 * VF->shade;			//	pn[ Index ] = 50 + 668 * cos( theta );
	SelectObject( bkhdc, br[ color ] );		// Index assigned brush

53366 shipped on 07 Jan 2019	Galago Pro
61396 shipped on 10 Dec 2019	Darter Pro
73783 shipped on 09 Dec 2020	Thelio
80429 shipped on 05 Jun 2021	Launch Keyboard

sudo chown root:root share  This really does make sense :)
https://www.3blue1brown.com/lessons/one-more-dim

#endif
#ifdef NOTES 
Plant genes:		->PL
	pollination, poison, spines, scent, Seed dispersal range:5 to 12??, camouflage?

Herbivore genes:	->HB
	Speed of motion, grazing range, scent, camouflage, eat spines, eat poison

Carnivore genes:	->CN
	speed of motion, attack range, camouflage?, sleep

	Separate speed from eat radius
moveCarnivore(), moveHerbivore() should use speed
eat() should use its own gene:: eat radius

Herbivores have a scent so HL->HB::scent vs carnivore CL->SN::smell  
#endif