//package ca;

import java.awt.Color;
import java.util.ArrayList;
import java.util.Random;

public class CellularAutomata
{       
	static Color[] colours = {Color.WHITE, Color.RED, Color.BLUE, Color.BLACK, Color.CYAN, Color.GREEN, Color.GRAY, Color.MAGENTA,
		Color.ORANGE, Color.PINK, Color.YELLOW};

	static Chart c;

	static Random randGen;
	static String spatial = "S";
	static String temporal = "bF";
	static int timeType = 2;
	static boolean metrics = false;

	static int totalPop = 0;
	static int[] population;
	static float[] bRate;
	static float[] dRate;
	static float[] bProb;
	static float[] dProb;	
	static int size;
	static int steps;
	static int maxTime;
	static float timestep;
	static double gillespTime = 0.0;

	static int[] grid;
	static int[][] neighbours;
	static int[] order;
	static double[][] cells;
	static double width;

	static ArrayList<int[]> changes; //NOT SURE THIS IS USED???

	static boolean[] countedCells;
	static boolean[] addedToPatchIndexes;
	static ArrayList<Integer> patchIndexes;
	
	static String rules;

	/**
	 * program should be started using the following settings:
	 * grid_length max_time time_step {inital_population birth_rate death_rate}*number of species.
	 * 
	 * optional arguments:
	 * 			-g  --> turn on the Grid 
	 * 			-c  --> turn on the Chart
	 * 			-m  --> turn on the spatial metrics (NP, LPI, Moran's I)
	 * 			-s x--> uses seed value x
	 * 			-r x--> runs with settings x where is can be either:
	 * 							"gil" for Gillespie simulator
	 * 						or
	 * 							update rule combination for CA in the form or Spatial,Temporal,TimeType from the options:
	 * 								Spatial:
	 * 									Sequential (S)
	 * 									Random (R)
	 * 								Temporal:
	 * 									birthFixed (bF)
	 * 									deathFixed(dF)
	 * 									birthConditional (bC)
	 * 									deathConditional (dC)
	 * 									singleRandom (sR)
	 * 									bothRandom (bR)
	 * 								Time Type:
	 * 									(1)
	 * 									(2)
	 */
	 
	public static void main(String[] args)
	{
	    long startTime = System.currentTimeMillis();
		int flags = 0;
		boolean GUI = false;
		boolean chart = false;
		boolean gillespie = false;
		long seed = System.currentTimeMillis();

		//determines the arguments of the program
		for(int i = 0; i < args.length; i++)
		{
			if(args[i].contains("-"))
			{
				flags++;
				if(args[i].equalsIgnoreCase("-g"))
				{
					GUI = true;
				}
				else if(args[i].equalsIgnoreCase("-c"))
				{
					chart = true;
				}
				else if (args[i].equalsIgnoreCase("-s"))
				{
					seed = Long.parseLong(args[i + 1]);
					flags++;
					i++;
				}
				else if(args[i].equalsIgnoreCase("-m"))
				{
					metrics = true;
				}
				else if (args[i].equalsIgnoreCase("-r"))
				{
					System.out.println(args[i+1]);
					rules = args[i+1];
					rules = rules.toLowerCase();

					i++;
					flags++;
				}
			}
		}

		// Do error checking on command line arguments
		if((args.length - 3 - flags) % 3 != 0)
		{
			System.out.println("Error with parameters. Please enter any flags " +
					"('-g' for no GUI, '-s' {long} for seed number, '-r' {rules} to change rules). " +
					"Required parameters are grid_length max_time time_step, and then for each species, the initial size, " +
					"the birth rate (>=0) and the death rate (>=0).");
			return;
		}
		if(((args.length - flags - 2)/3) > 10)
		{
			System.out.println("Too many species, only " + 10 + " allowed.");
			return;
		}

		// Set up some variables (why here - maybe because relies on knowing length of args ???)
		randGen = new Random(seed);
		population = new int[(args.length - flags - 2)/3];
		bRate = new float[(args.length - flags - 2)/3];
		dRate = new float[(args.length - flags - 2)/3];
		
		// Pick up command-line arguments and set up starting situation (sets totalPop, among other things)
		initialiseSpecies(args, flags);

		// Convert the rates to probabilities
		bProb = rate2probBin(bRate, timestep);
		dProb = rate2probBin(dRate, timestep);

		// Put populations onto the grid
		if(totalPop > (size*size))
		{
			System.out.println("There is not enough room on the grid for all of the species");
			return;
		}
		initiateGrid();
		initialiseNeighbours();

		// Set up some more variables ... mostly for measurements of output
		countedCells = new boolean[grid.length];
		addedToPatchIndexes = new boolean[grid.length];
		patchIndexes = new ArrayList<Integer>(grid.length);
		changes = new ArrayList<int[]>(grid.length);

		// Set up GUI and Chart
		if(GUI)
		{
			StdDraw.show(0);
			populateCells(grid);
			//done to make the GUI start up faster. Unsure why this stop the long initial load
			initiateGrid(new int[0]);
			initiateGrid(grid);
		}
		if(chart)
		{
			initiateChart();
		}

		// This section seems to set up an ordering to go through the grid. 
		// Obviously this is irrelevant if it turns out we've chosen to use the Gillespie ...
		order = new int[grid.length];

		//initiates the grid containing the spatial order
		for(int i = 0; i < order.length; i++)
		{
			order[i] = i;
		}
		if(spatial.equalsIgnoreCase("r"))
		{
			order = randomiseOrder(order);
		}
		else if(!(spatial.equalsIgnoreCase("s")))
		{
			System.out.println("Spatial input is incorrect. Must be either \"S\" or \"R\"");
			return;
		}

		// Step and time counters for all simulation rules
		double currTime = (double)0; // start at time zero
		steps = 0;
		// Output starting state at time 0
		printState(0);

		// -------------- GILLESPIE UPDATE LOOP ----------------
		
		if(rules.equalsIgnoreCase("gil"))
		{
			while( (currTime < maxTime) && (totalPop > 0) )			
			{
				// Reset local variables			
				double lambda = 0.0;		
				double[] birthRates = new double[bRate.length];
				double[] deathRates = new double[dRate.length];		
				double[] birthRatesNormalised = new double[bRate.length]; 
				double[] deathRatesNormalised = new double[dRate.length];
		
				double tau; // time until next event
				double eventProb; // for choosing our event
				double lowerBound = 0.0;
		
				// multiply by the population to get an overall rate.
				for(int i = 0; i < population.length; i++)
				{
					birthRates[i] = bRate[i] * population[i];
					deathRates[i] = dRate[i] * population[i];			
					lambda += birthRates[i] + deathRates[i];
				}		
				// Normalise the birth and death rates as a proportion of total rates
				for(int i = 0; i < birthRates.length; i++)
				{
					birthRatesNormalised[i] = birthRates[i]/lambda;
					deathRatesNormalised[i] = deathRates[i]/lambda;
				}
		
				// Determine time until next step
				tau = -Math.log(1 - randGen.nextDouble())/lambda; // time until next event
					
				// If tau is greater than timestep, then we need to print out intermed vals
				// and increment steps
				int tempCounter = 1;
				while (currTime + (tempCounter * timestep) < (int)(currTime + tau) ){
					printState((int)(currTime+tempCounter*timestep));
					steps++;
					tempCounter++;				
				}
				
				gillespTime += tau; // back into global variable ...
		
				// Choose next event
				eventProb = randGen.nextFloat();
				for(int i = 0; i < population.length; i++)
				{
					//checks all possible birth events
					if(eventProb >= lowerBound && eventProb < birthRatesNormalised[i] + lowerBound)
					{
						gillespieEvent(i, true);
						break;
					}
					lowerBound += birthRatesNormalised[i];
					//checks all possible death events
					if(eventProb >= lowerBound && eventProb < deathRatesNormalised[i] + lowerBound)
					{
						gillespieEvent(i, false);
						break;
					}
					lowerBound += deathRatesNormalised[i];
				
				}
				currTime = gillespTime; // pick this up from global variable affected by gillespieStep (!)
				
				// We've now done the updating, so output after new event here
				if (currTime > ((steps+1) * timestep)) {
					steps ++;
					printState((int)currTime);
				}
				
				// Update GUI and chart as required
				if(GUI)
				{
					updateGrid();
				}
				if(chart)
				{
					updateChart();
				}
				changes.clear();
			}			
		} else { // DOING A CA
		
		// CA - specific setup code
		int numSteps = (int)(maxTime / timestep);
		//creates variable which is pointer to the original grid for asynchronous updates
		int[] midstep = grid;
			
		// -------------- CA UPDATE LOOPS ----------------	
		
		// ----- SFB2
		if (rules.equalsIgnoreCase("sfb2")) {
				for (steps = 0; (steps < numSteps) && (totalPop > 0); steps++)
				{
					//bFixedStep(midstep);
					//midstep = births(grid, midstep); // do births and put into midstep
					midstep = birthsPoiss(grid, midstep);
					grid = deaths(grid, midstep); // do deaths onto midstep, scanning original grid (only originals can die)
					printState(steps+1);
					if(chart) updateChart();
					if(GUI) updateGrid();
					changes.clear();
				}
		} 
		else 
		{
			System.out.println("Incorrect string for rule. Only accepting sfb2 for the time being!");
			return;
		}	
		
		// FINAL OUTPUT --------
		//print statement for system going extinct
		if(totalPop == 0)
		{
			System.out.println("-System went extinct at step: " + steps);
		}
		System.out.println("-time: " + (System.currentTimeMillis() - startTime));
	} 
	
	}
	

	/************************* READ IN ARGUMENTS AND SET UP STARTING NUMBERS ***********************/
	 /**
	 * Obtains the initial values for the species, grid and number of steps.
	 * 
	 * @param args - arguments
	 * @param flags - no. of flags
	 */
	private static void initialiseSpecies(String[] args, int flags)
	{
		int counter = 0;

		String[] data = new String[args.length - flags];

		for(int i = 0; i < args.length; i++)
		{
			if(!args[i].contains("-"))
			{
				data[counter] = args[i];
				counter++;
			}
			else
			{
				if(args[i].equalsIgnoreCase("-s") || args[i].equalsIgnoreCase("-r"))
				{
					i++;
				}
			}
		}
		counter = 0;

		size = Integer.parseInt(data[0]);
		maxTime = Integer.parseInt(data[1]);
		timestep = Float.parseFloat(data[2]);

		for(int i = 3; i < data.length - 2; i++)
		{
			population[counter] = Integer.parseInt(data[i]);
			totalPop += population[counter];
			bRate[counter] = Float.parseFloat(data[i+1]);
			dRate[counter] = Float.parseFloat(data[i+2]);
			i += 2;
			counter++;
		}
	}

	/************************* POSITION ORIGINAL ORGANISMS ON GRID ***********************/
	/**
	 * Sets the initial locations for the cells of each species randomly
	 * 
	 * @param grid - grid to be placed
	 */
	private static void initiateGrid()
	{
		grid = new int[size*size];
		for(int i = 0; i < grid.length; i++)
		{
			grid[i] = 0;
		}

		int space = randGen.nextInt(grid.length);
		//for all species
		for(int i = 0; i < population.length; i++)
		{
			//for all organisms of that species
			for(int j = 0; j < population[i]; j++)
			{
				//until and empty cell is found
				while(grid[space] != 0)
				{
					space = randGen.nextInt(grid.length);
				}
				grid[space] = i + 1;
			}
		}
	}

	/************************* ??? ***********************/
	/**
	 * Sets the indexes for all of the neighbours for each site
	 */
	private static void initialiseNeighbours()
	{
		neighbours = new int[size * size][];

		for(int i = 0; i < neighbours.length; i++)
		{
			neighbours[i] = correctNumbers(i, neighbours.length);
		}
	}

	/************************* RANDOMISES ORDERING ?? ***********************/
	/**
	 * randomizes the order for the sites to be visited
	 * 
	 * @param order
	 * @return
	 */
	private static int[] randomiseOrder(int[] order)
	{
		int randNum;
		int temp;

		for(int i = 0; i < order.length; i++)
		{
			randNum = randGen.nextInt(order.length);
			temp = order[i];
			order[i] = order[randNum];
			order[randNum] = temp;
		}
		return order;
	}

	/************************* PRINT OUT CURRENT ITERATION ***********************/
	/**
	 * prints the current step, populations of each species and the total population of the system
	 *  also prints the LPI, NP and Moran's I is the -m flag is entered
	 * @param step
	 */
	private static void printState(int currTime)
	{
		System.out.print(currTime + " ");
		for(int i : population)
		{
			System.out.print(i + " ");
		}
		System.out.print(totalPop);
		if(metrics)
		{
			int[] patch = countPatches(grid);
			System.out.print(" " + patch[0] + " " + patch[1] +  " " + moransI(grid));
		}
		System.out.println(); //new line
		//printGrid(grid);
	}
	

	/************************* CHOOSES SITE AND EXECUTES BIRTH OR DEATH AT INDEX I - GILLESPIE ***********************/
	/**
	 * performs either a birth or death event on index i
	 * 
	 * @param i
	 * @param birth
	 */
	private static void gillespieEvent(int i, boolean birth)
	{
	    //System.out.println("gillespieEvent("+ i +","+ birth +") : "+ gillespTime);
		int cell = randGen.nextInt(grid.length);

		//keeps randomly selecting cells until one of the correct species is found
		while(grid[cell] == 0 || grid[cell] != (i + 1))
		{
			cell = randGen.nextInt(grid.length);
		}

		if(birth)
		{
			//performs a birth
			int[] hood = neighbours[cell];
			int birthSite = hood[randGen.nextInt(hood.length)];
			//System.out.print("species: "+ i +" cell: "+ cell +" birtSite : "+ birthSite +" > ");
			//for (int j=0;j<8;j++) System.out.print(hood[j] +" ");
			//System.out.println(" : "+ grid[birthSite]);
			if (grid[birthSite] == 0)
			{
				population[i]++;
				totalPop++;
				grid[birthSite] = i + 1;
				changes.add(new int[]{birthSite, i + 1});
			}
		}
		else
		{
			//performs a death
		    // System.out.println("species: "+ i +" cell: "+ cell +" kill ");
			grid[cell] = 0;
			population[i]--;
			totalPop--;
			changes.add(new int[] {cell, 0});
		}

	}

	/**************************************************************************/
	/************************* UPDATE RULES FOR CA HERE ***********************/
	/**************************************************************************/

	/************************* bFIXED - DO ONE ITERATION?? ***********************/
	/**
	 * performs the births for all of the cells and then the deaths
	 */
	 
	 // STARTING BY TRying to fix this one ...
	 
/*	private static void bFixedStep(int[] midstep)
	{
		if(timeType == 2)
		{
			//changes array to become its own grid if a synchronous option is chosen
			midstep = new int[grid.length];
		}
		midstep = births(grid, midstep);

		grid = deaths(grid, midstep);

	} */

	/**
	 * performs the deaths for all of the cells and then the births
	 */
	private static void dFixedStep(int[] midstep)
	{
		if(timeType == 2)
		{
			//changes array to become its own grid if a synchronous option is chosen
			midstep = new int[grid.length];
		}

		grid = deaths(grid, grid);

		grid = births(grid, midstep);
	}

	/**
	 * performs a birth, then death for each cell
	 */
	private static void bConditionalStep(int[] midstep)
	{
		if(timeType == 2)
		{
			//changes array to become its own grid if a synchronous option is chosen
			midstep = new int[grid.length];
		}

		//performs the checks in the specified order
		for(int i : order)
		{
			if(grid[i] != 0)
			{
				// birth step
				midstep[i] = grid[i];

				if(randGen.nextFloat() < bProb[grid[i] - 1])
				{
					int[] hood = neighbours[i];

					//randomly pick a neighoburing cell
					int birthSite = hood[randGen.nextInt(hood.length)];
					//birth only occurs if the chosen cell is empty in both the origninal and new grid
					if(grid[birthSite] == 0 && midstep[birthSite] == 0)
					{
						population[grid[i] - 1]++;
						totalPop++;
						midstep[birthSite] = grid[i];
						changes.add(new int[]{birthSite, grid[i]});
					}
				}

				//death step
				if(randGen.nextFloat() < dProb[grid[i] - 1])
				{
					population[grid[i] - 1]--;
					midstep[i] = 0;
					totalPop--;
					changes.add(new int[] {i, 0});
				}
			}
		}
		grid = midstep;
	}

	/**
	 * performs a death, then a birth for each cell
	 */
	private static void dConditionalStep(int[] midstep)
	{
		if(timeType == 2)
		{
			//changes array to become its own grid if a synchronous option is chosen
			midstep = new int[grid.length];
		}

		for(int i : order)
		{
			if(grid[i] != 0)
			{
				//death step

				if(randGen.nextFloat() < dProb[grid[i] - 1])
				{
					population[grid[i] - 1]--;
					midstep[i] = 0;
					totalPop--;
					changes.add(new int[] {i, 0});
				}
				else
				{
					// birth step

					midstep[i] = grid[i];                                   
					if(randGen.nextFloat() < bProb[grid[i] - 1])
					{
						int[] hood = neighbours[i];

						//randomly picks a neighbouring cell
						int birthSite = hood[randGen.nextInt(hood.length)];
						//only gives birth if the randomly selected cell is empty in both the original grid and the new one
						if(grid[birthSite] == 0 && midstep[birthSite] == 0)
						{
							population[grid[i] - 1]++;
							totalPop++;
							midstep[birthSite] = grid[i];
							changes.add(new int[]{birthSite, grid[i]});
						}
					}
				}
			}
		}
		grid = midstep;
	}

	/**
	 * performs a full step for a random setting, either both or single
	 * 
	 * @param both - true if bR, false if sR
	 * @param midstep
	 */
	private static void randomStep(boolean both, int[] midstep)
	{
		Boolean prev[] = new Boolean[order.length];
		Boolean second = false;

		if(timeType == 2)
		{
			//changes array to become its own grid if a synchronous option is chosen
			midstep = new int[grid.length];
		}

		for(int i = 0; i < prev.length; i++)
		{
			//an array which holds 1 if the event is randomly evaluated to be a birth, otherwise 0
			prev[i] = (randGen.nextInt(2) == 0);
		}

		for(int i = 0; i < order.length; i++)
		{
			int index = order[i];

			if(grid[index] != 0)
			{
				midstep[index] = grid[index];
				//if the event was earlier evalauted to be birth-first
				if(prev[i])
				{
					//performs a birth
					if(randGen.nextFloat() < bProb[grid[index] - 1])
					{
						int[] hood = neighbours[index];

						int birthSite = hood[randGen.nextInt(hood.length)];
						if(grid[birthSite] == 0 && midstep[birthSite] == 0)
						{
							population[grid[index] - 1]++;
							totalPop++;
							midstep[birthSite] = grid[index];
							changes.add(new int[]{birthSite, grid[index]});
						}
					}
					//if both events are to occur
					if(both)
					{
						//performs a death
						if(randGen.nextFloat() < dProb[grid[index] - 1])
						{
							population[grid[index] - 1]--;
							midstep[index] = 0;
							grid[index] = 0;
							totalPop--;
							changes.add(new int[] {index, 0});
						}
					}
					//if single
					else
					{
						//switches the value for the next iteration
						prev[i] = false;
					}
				}
				//if the event was earlier evaluated to be death-first
				else
				{
					//evaluated for death
					if(randGen.nextFloat() < dProb[grid[index] - 1])
					{
						population[grid[index] - 1]--;
						midstep[index] = 0;
						grid[index] = 0;
						totalPop--;
						changes.add(new int[] {index, 0});
					}
					//if death does not occur, evaulated for birth
					else
					{
						//if bR
						if(both)
						{
							if(randGen.nextFloat() < bProb[grid[index] - 1])
							{
								int[] hood = neighbours[index];

								int birthSite = hood[randGen.nextInt(hood.length)];
								if(grid[birthSite] == 0 && midstep[birthSite] == 0)
								{
									population[grid[index] - 1]++;
									totalPop++;
									midstep[birthSite] = grid[index];
									changes.add(new int[]{birthSite, grid[index]});
								}
							}
						}
						//if sR sets value to be true for next iteration
						else
						{
							prev[i] = true;
						}
					}
				}
			}
			//the grid has been fully iterated through for the first time, and sR was chosen (so needs to iterate twice)
			if((i == order.length - 1) && !second && !both)
			{
				//resets counter
				i = -1;
				second = true;
			}
		}
		//updates original grid for all changes
		grid = midstep;
	}

	/**
	 * performs the births for all of the cells in original grid, placing new data on the changed grid
	 * 
	 * @param orig
	 * @param changed
	 * @return
	 * I HAVE CHECKED THIS AND THINK IT'S OK, at least for sfb2
	 * Scans orig (A) and puts all of these organisms into changed, plus any births where site empty in A and B
	 */
	private static int[] births(int[] orig, int[] changed)
	{
		//goes in the defined order
		for(int i: order)
		{
			if(orig[i] != 0) // if organism at site
			{
				//this cell survives this step no matter what
				changed[i] = orig[i];

				if((randGen.nextFloat() < bProb[orig[i] - 1]))
				{
					int[] hood = neighbours[i];
					//picks random neighbour
					int birthSite = hood[randGen.nextInt(hood.length)];
					// only gives birth if the selected neighbour is free in both the original and new grid
					// this means that being early in the sequence is advantageous
					if(orig[birthSite] == 0 && changed[birthSite] == 0)
					{
						population[orig[i]- 1]++;
						totalPop++;
						changed[birthSite] = orig[i];
						changes.add(new int[]{birthSite, orig[i]});
					}
				}

			}
		}
		//System.out.println(changed.length);
		return changed;
	}
	
	/**
	 * Performs the births for all of the cells in original grid, placing new data on the changed grid
	 * Allows multiple births for each organism by converting rates to probabilities using Poisson dist
	 *
	 * @param orig
	 * @param changed
	 * @return
	 * I HAVE CHECKED THIS AND THINK IT'S OK, at least for sfb2
	 * Scans orig (A) and puts all of these organisms into changed, plus any births where site empty in A and B
	 */
	private static int[] birthsPoiss(int[] orig, int[] changed)
	{
		int numBirths;
		//goes in the defined order
		for(int i: order)
		{
			if(orig[i] != 0) // if organism at site
			{
				//this cell survives this step no matter what
				changed[i] = orig[i];
				
				// work out number of births
				// OK System.out.println("Birth rate at current site " + bRate[orig[i]-1]);
				numBirths = numberOfBirths(bRate[orig[i]-1] * timestep);
				System.out.println("Expecting to generate " + numBirths + " births.");
				
				if((randGen.nextFloat() < bProb[orig[i] - 1]))
				{
					int[] hood = neighbours[i];
					//picks random neighbour
					int birthSite = hood[randGen.nextInt(hood.length)];
					// only gives birth if the selected neighbour is free in both the original and new grid
					// this means that being early in the sequence is advantageous
					if(orig[birthSite] == 0 && changed[birthSite] == 0)
					{
						population[orig[i]- 1]++;
						totalPop++;
						changed[birthSite] = orig[i];
						changes.add(new int[]{birthSite, orig[i]});
					}
				}

			}
		}
		//System.out.println(changed.length);
		return changed;
	}

	/**
	 * performs the deaths for all of the cells in original grid, placing new data on the changed grid
	 * 
	 * @param orig
	 * @param changed
	 * @return
	 * Scans grid (A) and makes changes in changed (B), returning changed
	 */
	private static int[] deaths(int[] grid, int[] changed)
	{
		//performs the checks in the specified order
		for(int i: order)
		{
			//System.out.println(i);
			if(grid[i] != 0)
			{
				if((randGen.nextFloat() < dProb[grid[i] - 1]))
				{
					population[grid[i] - 1]--;
					changed[i] = 0;
					totalPop--;
					changes.add(new int[] {i, 0});
				}
			}
		}

		return changed;
	}

	/**************************************************************************/
	/************************* END OF UPDATE RULES ****************************/
	/**************************************************************************/

	/************************* GET NEIGHBOURS OF CELL ***********************/	
	/**
	 * returns a list of all of the correct indexes for the neighbours of the provided cell (i)
	 * 
	 * @param i - index of the cell to find neighbours for
	 * @param length - size of the array
	 * @returns a list of indexes for the grid for the neighbours of the given cell
	 */
	private static int[] correctNumbers(int i, int length)
	{
		int[] neighbourIndex = new int[8];
		int rowLength = (int)Math.sqrt(length);

		// for the top-left neighbour
		neighbourIndex[0] = (i - rowLength - 1);
		if(i < rowLength)
		{
			if(i == 0)
			{
				neighbourIndex[0] = length - 1; 
			}
			else
			{
				neighbourIndex[0] += length;
			}
		}
		else if(i % rowLength == 0)
		{
			neighbourIndex[0] = i - 1;
		}


		// for the neighbour directly above 
		neighbourIndex[1] = (i - rowLength);
		if(i < rowLength)
		{
			neighbourIndex[1] += length;
		}

		//for the top-right neighbour
		neighbourIndex[2] = (i - rowLength + 1);
		if(i < rowLength)
		{
			if(i == rowLength - 1)
			{
				neighbourIndex[2] = length - rowLength;
			}
			else
			{
				neighbourIndex[2] += length;
			}
		}
		else if((i + 1) % rowLength == 0)
		{
			neighbourIndex[2] -= rowLength;
		}

		// for the left neighbour
		neighbourIndex[3] = i - 1;
		if(i % rowLength == 0)
		{
			neighbourIndex[3] += rowLength;
		}

		// for the right neighbour
		neighbourIndex[4] = i + 1;
		if((i + 1) % rowLength == 0)
		{
			neighbourIndex[4] -= rowLength;
		}

		// for the bottom-left neighbour
		neighbourIndex[5] = (i + rowLength - 1);
		if(i > length - rowLength - 1)
		{
			if(i == (length - rowLength))
			{
				neighbourIndex[5] = rowLength - 1;
			}
			else
			{
				neighbourIndex[5] -= length;
			}
		}
		else if(i % rowLength == 0)
		{
			neighbourIndex[5] += rowLength; 
		}

		// for the neighbour directly bellow
		neighbourIndex[6] = i + rowLength;
		if(i > length - rowLength - 1)
		{
			neighbourIndex[6] -= length;
		}

		// for the bottom-right neighbour
		neighbourIndex[7] = i + rowLength + 1;
		if(i > length - rowLength - 1)
		{
			if(i == length - 1)
			{
				neighbourIndex[7] = 0;
			}
			else
			{
				neighbourIndex[7] -= length;
			}
		}
		else if((i + 1) % rowLength == 0)
		{
			neighbourIndex[7] -= rowLength;
		}

		return neighbourIndex;
	}

	/************************* PRINT GRID TO STANDARD OUTPUT ***********************/
	/**
	 * textual print-out of the grid
	 * 
	 * @param grid
	 */
	private static void printGrid(int[] grid)
	{
		int l = (int) Math.sqrt(grid.length);
		for(int i = 0; i < grid.length; i++)
		{

			System.out.print(grid[i] + " ");

			if((i + 1)%l == 0)
			{
				System.out.println("");
			}
		}
		System.out.println("");
	}
	
	/************************* MORAN'S I FOR GIVEN GRID ***********************/
	/**
	 * calculates moran's I for the given grid
	 * @param grid
	 * @return
	 */
	private static float moransI(int[] grid)
	{
		float moransI = 0;
		float xBar = 0;
		float sac = 0;
		float tv = 0;

		for(int i:population)
		{
			xBar += i;
		}
		xBar /= grid.length;

		//for all cells
		for(int i = 0; i < grid.length; i++)
		{
			//for each of its neighbours
			for(int j:neighbours[i])
			{
				if(grid[i] == grid[j])
				{
					//for two empty cells
					if(grid[i] == 0)
					{
						sac += (0-xBar)*(0-xBar);
					}
					//for two cells of the same species
					else
					{
						sac += (1-xBar)*(1-xBar);
					}
				}
				//one empty cell and a non-matching cell
				else
				{
					sac += (0-xBar)*(1-xBar);
				}
			}
			if(grid[i] == 0)
			{
				tv += (0 - xBar)*(0 - xBar);
			}
			else
			{
				tv += (1 - xBar)*(1 - xBar);
			}
		}

		float n = grid.length;
		float connections = n*neighbours[0].length;

		moransI = (n/connections) * (sac/tv);

		return moransI;
	}

	/************************* NUMBER PATCHES AND LARGEST PATCH INDEX ***********************/
	/**
	 * counts the number of patches and the largest patch index
	 * 
	 * @param grid
	 * @returns the number of patches, largest patch size
	 */
	private static int[] countPatches(int[] grid)
	{
		int patches = 0;
		int biggestPatch = 0, currentPatch = 0;

		//initiates the two boolean arrays
		for(int i = 0; i < countedCells.length; i++)
		{
			countedCells[i] = false;
			addedToPatchIndexes[i] = false;
		}

		for(int i = 0; i < grid.length; i++)
		{
			//all non-empty, uncounted cells
			if(!countedCells[i] && grid[i] != 0)
			{
				currentPatch = 0;
				//add to arrayList
				patchIndexes.add(i);
				addedToPatchIndexes[i] = true;
				//for all neighbours
				while(!patchIndexes.isEmpty())
				{
					int currentIndex = patchIndexes.remove(0);
					currentPatch++;
					countedCells[currentIndex] = true;
					addedToPatchIndexes[currentIndex] = false;
					for(int n:neighbours[currentIndex])
					{
						if(grid[n] != 0 && grid[n] == grid[i] && !countedCells[n] && !addedToPatchIndexes[n])
						{
							patchIndexes.add(n);
							addedToPatchIndexes[n] = true;
						}
					}
				}

				patches++;
				//checks for new biggest
				if (currentPatch > biggestPatch)
				{
					biggestPatch = currentPatch;
				}

			}
		}

		return new int[]{patches, biggestPatch};
	}

	/************************* ??????????? ***********************/
	/**
	 * populates the array cells with the coordinates for each cell
	 * 
	 * @param grid
	 */
	public static void populateCells(int[] grid)
	{
		int rowLength = (int)Math.sqrt(grid.length);
		int i = 0;
		double incr = ((double)1/rowLength) * (StdDraw.xmax - StdDraw.xmin);
		width = incr/2;
		cells = new double[grid.length][];

		double y = 0;
		double x = 0;

		for(int row = 0; row < (int)Math.sqrt(grid.length); row++)
		{
			for(int column = 0; column < (int)Math.sqrt(grid.length); column++)
			{
				cells[i] = new double[]{x, y};
				i++;
				x += incr;
			}
			y += incr;
			x = 0;
		}
		StdDraw.show(0);
	}

	/************************* OUTPUT CODE ***********************/

	public static void initiateChart()
	{
		c = new Chart(population);;
	}

	public static void updateChart()
	{
		c.updatePopulations(population);
	}

	/**
	 * displays the initial layout of the system
	 * 
	 * @param grid
	 */
	public static void initiateGrid(int[] grid)
	{
		for(int i = 0; i < grid.length; i++)
		{
			StdDraw.setPenColor(colours[grid[i]]);
			StdDraw.filledSquare(cells[i][0], cells[i][1], width);
		}
	}

	/**
	 * updates the display to show which cells have changed 
	 */
	public static void updateGrid()
	{
		while(changes.size() != 0)
		{
			int[] change = changes.remove(0);

			StdDraw.setPenColor(colours[(change[1])]);
			if(StdDraw.getPenColor() == Color.WHITE)
			{
				StdDraw.filledSquare(cells[change[0]][0], cells[change[0]][1], width + 0.001);
			}
			else
			{
				StdDraw.filledSquare(cells[change[0]][0], cells[change[0]][1], width);
			}
		}
		StdDraw.show(1);
	}
	
	/************************ CODE TO CONVERT RATES TO PROBABILITIES *************************/

	public static float[] rate2probBin(float[] rate, double timeunit) { 
		float[] proba;
		proba = new float[rate.length];
		
		for (int i=0; i < rate.length; i++) 
		{ 
			proba[i] = ((float)1 - (float)Math.exp(-rate[i]*timeunit));
       	}
       	return proba;		
    }
    
    /* Draws from a Poisson distribution a number of births to execute on the basis of a rate and timeunit */
    /*public static int rate2numPoiss(float[] rate, double timeunit) { 
		int num = 0;
		float lambda = rate*timeunit;
		//num = Poisson.random(rate*timeunit);
       	RandomEngine engine = new DRand();
		Poisson poisson = new Poisson(lambda, engine);
		int poissonObs = poisson.nextInt();
       	return num;		
    }*/
    
    /* This is Pat's function from CAR in ca7.1*/
    public static int numberOfBirths(double lambda){
	double L     = 0.0;
	double U     = randGen.nextDouble();
//	double U     = gen.nextDouble();
	double factB = 1.0;
	for (int b=0;b<8;b++){
	    if (b > 0) factB = factB * b;
	    L = L + Math.pow(lambda*timestep,(double)b) * Math.pow(Math.E,-lambda*timestep)/factB;
	    if (L >= U) return b;
	}
	return 8; // size of Moore neighbourhood
    }
    
    /* Generates a matrix where size of slice corresponding to each ???????
    public static float[][] poissDist(float[] rate, double timeunit) {
    	float[] proba;
    	proba = new float[rate.length][8]; // first dimension for rates for each species, second for each neighbour
    	
    	float lambda = rate*timeunit; // parameters of poiss distribution
    	 
    	for (int i=0; i < rate.length; i++)
    	{
    		for (int j=0; j < 8; j++)
  			{
  				// stuff here
  				proba[i][j] = ( Math.pow(lambda[i],j)  / factorial(j) ) * Math.exp(-lambda[i]);
  			}
    	}
    } */
    
    public static int factorial(int n) {
        int fact = 1; // this  will be the result
        for (int i = 1; i <= n; i++) {
            fact *= i;
        }
        return fact;
    }
   
}

/*		//determines the exact setting type before looping to reduce expensive "if" statements for every step
		else
		{
			

			if(tempor.equalsIgnoreCase("f"))
			{
				if(ord.equalsIgnoreCase("b"))
				{
					//System.out.println("steps=" + steps);
					//for(; steps < i && totalPop > 0; steps++)
					for (steps = 0; (steps < numSteps) && (totalPop > 0); steps++)
					{
						bFixedStep(midstep);
						printState(steps+1);
						if(GUI)
						{
							updateGrid();
						}
						if(chart)
						{
							updateChart();
						}
						changes.clear();
						if(spatial.equalsIgnoreCase("r"))
						{
							//re-randomises the spatial order
							order = randomiseOrder(order);
						}
					}
				}
				else if(ord.equalsIgnoreCase("d"))
				{
					for(; steps < i && totalPop > 0; steps++)
					{
						dFixedStep(midstep);
						printState(steps);
						if(GUI)
						{
							updateGrid();
						}
						if(chart)
						{
							updateChart();
						}
						changes.clear();
						if(spatial.equalsIgnoreCase("r"))
						{
							order = randomiseOrder(order);
						}
					}
				}
				else
				{
					temporError = true;
				}
			}
			else if(tempor.equalsIgnoreCase("c"))
			{
				if(ord.equalsIgnoreCase("b"))
				{
					for(; steps < i && totalPop > 0; steps++)
					{
						bConditionalStep(midstep);
						printState(steps);
						if(GUI)
						{
							updateGrid();
						}
						if(chart)
						{
							updateChart();
						}
						changes.clear();
						if(spatial.equalsIgnoreCase("r"))
						{
							order = randomiseOrder(order);
						}
					}
				}
				else if(ord.equalsIgnoreCase("d"))
				{
					for(; steps < i && totalPop > 0; steps++)
					{
						dConditionalStep(midstep);
						printState(steps);
						if(GUI)
						{
							updateGrid();
						}
						if(chart)
						{
							updateChart();
						}
						changes.clear();
						if(spatial.equalsIgnoreCase("r"))
						{
							order = randomiseOrder(order);
						}
					}
				}
				else
				{
					temporError = true;
				}
			}
			else if(tempor.equalsIgnoreCase("r"))
			{
				if(ord.equalsIgnoreCase("b"))
				{
					for(; steps < i && totalPop > 0; steps++)
					{
						randomStep(true, midstep);
						printState(steps);
						if(GUI)
						{
							updateGrid();
						}
						if(chart)
						{
							updateChart();
						}
						changes.clear();
						if(spatial.equalsIgnoreCase("r"))
						{
							order = randomiseOrder(order);
						}
					}
				}
				else if(ord.equalsIgnoreCase("s"))
				{
					for(; steps < i && totalPop > 0; steps++)
					{
						randomStep(false, midstep);
						printState(steps);
						if(GUI)
						{
							updateGrid();
						}
						if(chart)
						{
							updateChart();
						}
						changes.clear();
						if(spatial.equalsIgnoreCase("r"))
						{
							order = randomiseOrder(order);
						}
					}
				}
				else
				{
					temporError = true;
				}
			}
			else
			{
				temporError = true;
			}

			if(temporError)
			{
				System.out.println("Temporal input is incorrect. Must be either \"bF\", \"dF\", \"bC\", \"dC\", " +
						"\bR\" or \"sR\"");
				return;
			}
		}
*/

	/**
	 * prints the current step, populations of each species and the total population of the system
	 *  also prints the LPI, NP and Moran's I is the -m flag is entered
	 * @param step
	 */ 
	 // OLD VERSION ...
	/*private static void printState(int step)
	{
		System.out.print(step + 1 + " ");
		for(int i : population)
		{
			System.out.print(i + " ");
		}
		if(metrics)
		{
			int[] patch = countPatches(grid);
			System.out.println(totalPop + " | " + patch[0] + " " + patch[1] +  " " + moransI(grid));
		}
		else
		{
			System.out.println(totalPop);
		}
		//printGrid(grid);
	} */

	/************************* EXECUTES ONE GILLESPIE EVENT ***********************/
	/**
	 * performs each of the gillespie steps
	 */
	/*private static void gillespieStep()
	{
		double lambda = 0.0;

		double[] birthRates = new double[bRate.length];
		double[] deathRates = new double[dRate.length];

		double[] birthRatesNormalised = new double[bRate.length]; 
		double[] deathRatesNormalised = new double[dRate.length];

		double tau;
		double eventProb; // for choosing our event
		double lowerBound = 0.0;

		for(int i = 0; i < population.length; i++)
		{
			//math function figures out the average rate from the probability using a negative exponential distribution
			//then multiplies by the population to get an overall rate.
			birthRates[i] = bRate[i] * population[i];
			deathRates[i] = dRate[i] * population[i];			
			lambda += birthRates[i] + deathRates[i];
		}

		// Normalise the birth and death rates as a proportion of total rates
		for(int i = 0; i < birthRates.length; i++)
		{
			birthRatesNormalised[i] = birthRates[i]/lambda;
			deathRatesNormalised[i] = deathRates[i]/lambda;
		}

		// Determine time until next step
		tau = -Math.log(1 - randGen.nextDouble())/lambda; // time until next event
		gillespTime += tau; // back into global variable ...
	
		// Choose next event
		eventProb = randGen.nextFloat();
		for(int i = 0; i < population.length; i++)
		{
			//checks all possible birth events
			if(eventProb >= lowerBound && eventProb < birthRatesNormalised[i] + lowerBound)
			{
				gillespieEvent(i, true);
				break;
			}
			lowerBound += birthRatesNormalised[i];
			//checks all possible death events
			if(eventProb >= lowerBound && eventProb < deathRatesNormalised[i] + lowerBound)
			{
				gillespieEvent(i, false);
				break;
			}
			lowerBound += deathRatesNormalised[i];
		}
	}
	*/

					/*System.out.println(rules);
					if(rules.equalsIgnoreCase("gil"))
					{
						gillespie = true;
					}
					else
					{
						spatial = rules.substring(0, 1);
						temporal = rules.substring(1, 3);
						timeType = Integer.parseInt(rules.substring(3));
					}*/
					
					
					/*		String tempor = temporal.substring(temporal.length() - 1);
		String ord = temporal.substring(0, 1);
		boolean temporError = false; */