## set user-defined parameters

# arena size, should be 8<=s<=60, best at about 50
s = 50

# birth probability, must be 0<b<1, best at 0.1 to 0.8
b = 1/3

# death probability, must be 0<=d<b, best at 0 to b/4
d = 0

# number of trials, must be 1, 2, 3, or 4, larger means more data and more runtime
N = 4

# starting setup: 1 for center and 2 for corners
startup = 1

# curve fit: 0 for no, 1 for fitting r with theoretical K, 2 for 2-parameter fit
curvefit = 2

## set fixed parameters

n = s
T = 100
K = n^2
plotz = 2

## initialize multirun data structures

x = c()
y = c()

## start a run

for (run in 1:N)
{

## define and initialize data structures

# n rows are i=2 to i=n+1, similarly for columns
M <- matrix(0,nrow=n+2,ncol=n+2)

P <- rep(0,T+1)


# if startup==1 mark central squares as occupied (M=5)
if (startup==1)
	for (i in 1:2)
		for (j in 1:2)
			M[n/2+i,n/2+j] <- 5

# if startup==2, mark corners as occupied (M=5)
if (startup==2)
	for (i in c(2,n+1))
		for (j in c(2,n+1))
			M[i,j] <- 5


P[1] = 4

## run first step -- no deaths

# mark as available (M=s) if adjacent to s occupied
for (i in 2:(n+1))
	for (j in 2:(n+1))
		if (M[i,j]<5)
			{
			s = 0
			if (M[i-1,j]==5)
				s = s+1
			if (M[i+1,j]==5)
				s = s+1
			if (M[i,j-1]==5)
				s = s+1
			if (M[i,j+1]==5)
				s = s+1
			M[i,j] = s
			}

# occupy available sites with probability b for each occupied neighbor
for (i in 2:(n+1))
	for (j in 2:(n+1))
		if (M[i,j]<5)
			if (M[i,j]>0)
				if (min(runif(M[i,j]))<b)
					M[i,j] = 5

P[2] = length(which(M==5))

## run simulation

for (t in 2:T)
	{
# unoccupy filled sites with probability d
	for (i in 2:(n+1))
		for (j in 2:(n+1))
			if (M[i,j]==5)
				if (runif(1)<d)
					M[i,j] = 0

# mark as available (M=s) if adjacent to s occupied
	for (i in 2:(n+1))
		for (j in 2:(n+1))
			if (M[i,j]<5)
				{
				s = 0
				if (M[i-1,j]==5)
					s = s+1
				if (M[i+1,j]==5)
					s = s+1
				if (M[i,j-1]==5)
					s = s+1
				if (M[i,j+1]==5)
					s = s+1
				M[i,j] = s
				}

# occupy available sites with probability b for each occupied neighbor
	for (i in 2:(n+1))
		for (j in 2:(n+1))
			if (M[i,j]<5)
				if (M[i,j]>0)
					if (min(runif(M[i,j]))<b)
						M[i,j] = 5

	P[t+1] = length(which(M==5))
	}

## collect results

x1 = P
T = min(T,which(P>=K))
T1 = T
dP = P[2:(T+1)]-P[1:T]
y1 = dP
T = min(T,which(dP<0))
T2 = T

x = c(x,P[1:T])
y = c(y,dP[1:T])

## end the run loop

}


## define functions for data analysis

# Define the function f(p,x) for the model y=qf(p,x), where p is the nonlinear parameter.
f <- function(p,x) x*(p-x)

# Define the function F(p) that determines RSS, given optimal q.
F <- function(p) -(sum(y*f(p,x)))^2/sum((f(p,x))^2)	


## analyze the data

numpoints = length(x)

# find 1-parameter best fit

X = f(K,x)
r1 = sum(X*y)/sum(X^2)
RSS1 = sum(y^2)+F(K)
AIC1 = numpoints*log(RSS1/numpoints)+2*2

# find 2-parameter best fit

result <- optimize(f=F, lower=0.5*K, upper=2*K, maximum=FALSE)
K2 <- result$minimum
X = f(K2,x)
r2 <- sum(X*y)/sum(X^2)
RSS2 <- sum(y^2)+result$objective
AIC2 = numpoints*log(RSS2/numpoints)+2*3

# plot data and models

par(mfrow=c(2,2))

#1

plot(0:T,P[1:(T+1)],xlab="time",ylab="population",type="l")

#2

plot(1:T,dP[1:T],xlab="time",ylab="change in population")

#3

plot(x,y,xlab="population",ylab="change in population")

X = 0:K2

Y1 = r1*f(K,X)
if (curvefit==1)
	lines(X,Y1,col="blue")

Y2 = r2*f(K2,X)
if (curvefit==2)
	lines(X,Y2,col="red")

#4

i4 = which(x>6)
plot(x,y/x,xlab="population",ylab="relative change (dP/P)",ylim=c(0,2*r2*K2))

X[1] = 1
Y1[1] = r1*K
Y2[1] = r2*K2

if (curvefit==1)
	lines(X,Y1/X,col="blue")

Y2 = r2*f(K2,X)
if (curvefit==2)
	lines(X,Y2/X,col="red")

# report results for chosen fitting protocol

K
r1
K2
r2