Model Thinking- Game of Life- varying starting densities

GameOfLifeStartingDensities.r

# this code does a little experiment to see how the proportion living in a 500x500 game of life board changes 
# over time by starting population density.

# we first define 2 functions (from Vadim Vinichenko) for the game of life. 
# His functions are convenient because they scale well.

shiftMatrix <- function(mx, dr, dc) {
	#Shift the matrix by dr (delta r) rows and dc columns
	#by adding e.g. dr rows of zeros and removing dr rows from the other side 
	
	nr <- nrow(mx)
	nc <- ncol(mx)
	
	#If the matrix is shifted by more than its nrow or ncol, we get a matrix of zeros 
	if (abs(dr) >= nr || abs(dc) >= nc) {
		mx <- matrix(0, nrow = nr, ncol = nc)
		return(mx)
	}
	
	#Rows:
	if (dr > 0) {
		mx <- rbind(mat.or.vec(dr, nc), mx)
		mx <- mx[1:nr,]
	} else if (dr < 0) {
		mx <- rbind(mx, mat.or.vec(-dr, nc))
		mx <- mx[(1 - dr):(nr - dr),]
	}
	
	#Columns:
	if (dc > 0) {
		mx <- cbind(mat.or.vec(nr, dc), mx)
		mx <- mx[,1:nc]
	} else if (dc < 0) {
		mx <- cbind(mx, mat.or.vec(nr, -dc))
		mx <- mx[,(1 - dc):(nc - dc)]
	}
	return(mx)
}

life_cycle <- function(mx) {
	#Move the board one generation forward
	
	mx0 <- matrix(0, nrow = nrow(mx), ncol = ncol(mx))
	
	#Produce 8 "shifted" boards and add them up
	for (n in (-1:1)) {
		for (m in (-1:1)) {
			if (n !=0 || m !=0) {
				mx0 <- mx0 + shiftMatrix(mx, n, m)
			}
		}
	}
	
	#Deaths and births
	mx[mx0 > 3 | mx0 < 2] <- 0
	mx[mx0 == 3] <- 1
	return(mx)
}

# begin experiment:

# set parameters:
n_rows <- 500
n_cols <- 500
n_cells <- n_rows * n_cols
n_cycles <- 100
sleep_time <- 0.1
clrs <- c("#f0f0f0", "#2f81c1") #colors for empty squares and cells
rnd_threshold <- seq(.01,.5,by=.005) # 0 - empty board; 1 - all squares are filled

# this will be the container for results:
prop_on <- matrix(ncol=length(rnd_threshold),nrow=(n_cycles+1))

# the test: (takes a long time to run! ca 15-30 minutes)
for (z in 1:length(rnd_threshold)){
	# each iteration starts with the same board
	set.seed(1)
	board <- matrix(0, nrow = n_rows, ncol = n_cols)
	# but the density of 'on' cells changes 
	board[runif(n_cells, 0, 1) < rnd_threshold[z]] <- 1
	prop_on[1,z] <- sum(board)/n_cells
	# here we iterate n_cycles times for the given density parameter
	for (i in (1:n_cycles)) {
		board <- life_cycle(board)
		prop_on[(i+1),z] <- sum(board)/n_cells
	}
}

# save results for later analysis if desired:
#save(prop_on,file="proportionOnStartingDensityTests.Rdata")
#load("proportionOnStartingDensityTests.Rdata")

# look at results:

# first graphic requires "fields" package:
# install.packages("fields")
#png("GameOfLifeSurface.png")
fields:::image.plot(x=0:n_cycles,y=rnd_threshold,prop_on,xlab="time",ylab="starting population density density",
		main="population density over time by starting population density")
#dev.off()

# second graphic:
#png("GameOfLifeRatioLivingByStartingDensity.png")
plot(x=rnd_threshold,y=prop_on[101,]/prop_on[1,],
		type='l',
		main="proportion of original livings alive after 100 iterations\nby starting density",
		ylab="ratio of t=100 density to t=0 density",
		xlab="starting density")
#dev.off()

# third graphic:
#png("GameOfLifeProportionLiving100IterationsbyStartingDensity.png")
plot(x=rnd_threshold,y=prop_on[101,],
		type='l',
		main="proportion living after 100 iterations by starting density",
		xlab="starting density",
		ylab="proportion of total cells living at t=100")
#dev.off()

# end