Model Thinking- Aggregation- Game of Life- Comparing dead and wrapped borders by starting density

GOL_startingdensbyborderassumption500x500_100iterations.r

# FIXED!!

# 2 simulations run, both take a long time. This code will dish it out to 8 cores.
# I will post the 2 data objects produced as .Rdata files somewhere so you can just load them
# and follow the code to reproduce the figures, etc.

# rnd_threshold are starting densities:
rnd_threshold <- seq(.01,.99,by=.01)

ParThreshold <- function(thr,n_rows=500,n_cols=500,n_cycles=5){
	
	# the 2 GOL functions from Vadim Vinichenko:
	# recall the assumption here that borders are dead, rather than wrapping
	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)
	}
	
	lifeCycle <- 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)
	}
	
	# produce a starting board:
	n_cells <- n_rows * n_cols
	# assures same random number sequence
	set.seed(1)
	board <- matrix(0, nrow = n_rows, ncol = n_cols)
	# but the density of 'on' cells changes with 'thr' (threshold)
	board[runif(n_cells, 0, 1) < thr] <- 1
	# prop_on is the output (proportion living after each iteration
	prop_on <- vector(length = n_cycles+1)
	prop_on[1] <- sum(board)/n_cells
	# here we iterate n_cycles times for the given density parameter
	for (i in (1:n_cycles)) {
		board <- lifeCycle(board)
		prop_on[(i+1)] <- sum(board)/n_cells
	}
	return(prop_on)
}

# This is a big task, so we'll dish it out to 8 available server cores
# really THINK TWICE before running this code. Anyway it'll only work if you have a 
# 8+ core machine at your disposal. You could run it on a single core using sapply()
# but you'd be waiting possibly more than a couple days for it to finish...
library(parallel)
cl <- makeCluster(getOption("cl.cores", 8))
A <- parSapply(cl,rnd_threshold,ParThreshold,n_cycles=1000)
stopCluster(cl)

# and this little line saves the output in a file that I will attach
# and that all can analyze as they see fit..
colnames(A) <- rnd_threshold
rownames(A) <- 0:1000
save(A,file="H:/compartit/triffe/GOL_Dens1000Dead.Rdata")

fields:::image.plot(y=rnd_threshold,x=0:149,A[1:150,],main="Game of Life, 500x500\nInitial Density Test with dead borders",xlab="time",ylab="density")
# note that the tail isn't dead:
fields:::image.plot(A[900:1001,])
# -----------------------
# make same object using border wrapping function instead
ParThreshold2 <- function(thr,n_rows=500,n_cols=500,n_cycles=5){
	
	GOLborderwrap <- function(M){
		nc <- ncol(M)
		nr <- nrow(M)
		# augmented M wraps around edges to form a taurus shape-
		# i.e. a continuous world rather than dead borders.
		M.aug <- cbind(c(M[nr,nc],M[,nc],M[1,nc]),
				rbind(M[nr,],M,M[1,]),
				c(M[nr,1],M[,1],M[1,1]))
		
		# how many neighbors?:
		# left side
		M.out <-	M.aug[1:nr, 1:nc] +
				M.aug[2:(nr+1), 1:nc] +
				M.aug[3:(nr+2), 1:nc] +
				
				# middle:
				M.aug[1:nr,2:(nc+1)] +
				M.aug[3:(nr+2),2:(nc+1)] +
				
				# right:
				M.aug[1:nr,3:(nc+2)] +
				M.aug[2:(nr+1),3:(nc+2)] +
				M.aug[3:(nr+2),3:(nc+2)] 
		
		# birth and death rules:
		M[M.out > 3 | M.out < 2] <- 0
		M[M.out == 3 | M.out == 2] <- 1
		return(M)
	}
	
	# produce a starting board:
	n_cells <- n_rows * n_cols
	# assures same random number sequence
	set.seed(1)
	board <- matrix(0, nrow = n_rows, ncol = n_cols)
	# but the density of 'on' cells changes with 'thr' (threshold)
	board[runif(n_cells, 0, 1) < thr] <- 1
	# prop_on is the output (proportion living after each iteration
	prop_on <- vector(length = n_cycles+1)
	prop_on[1] <- sum(board)/n_cells
	# here we iterate n_cycles times for the given density parameter
	for (i in (1:n_cycles)) {
		board <- GOLborderwrap(board)
		prop_on[(i+1)] <- sum(board)/n_cells
	}
	return(prop_on)
}

cl <- makeCluster(getOption("cl.cores", 8))
B <- parSapply(cl,rnd_threshold,ParThreshold2,n_cycles=1000)
stopCluster(cl)
colnames(B) <- rnd_threshold
rownames(B) <- 0:1000
save(B,file="H:/compartit/triffe/GOL_Dens1000Wrap.Rdata")
#-------------------------------------------------------
#plot(0:1000,B[,1],type='l',ylim=c(0,1),col="#80808050")
#for (i in 2:99){
#	lines(0:1000,B[,i],col="#80808050")
#}

# you can skip the above and just load these 2 objects to continue.
# you can download the above 2 objects (A and B) here:
https://sites.google.com/site/timriffepersonal/r-code/GOL_Dens1000Dead.Rdata?attredirects=0&d=1
https://sites.google.com/site/timriffepersonal/r-code/GOL_Dens1000Wrap.Rdata?attredirects=0&d=1
# ...Dead is A, and ...Wrap is B

# load in the data objects:
load("H:/compartit/triffe/GOL_Dens1000Dead.Rdata")
load("H:/compartit/triffe/GOL_Dens1000Wrap.Rdata")

# http://s19.postimage.org/n63hhdh2r/GOLdensity_Tests_Deadvs_Wrapping500board.png
# comapre densities over first 150 iterations (surfaces):
# install.packages("fields") # if necessary
par(mfrow=c(1,2))
fields:::image.plot(y=rnd_threshold,x=0:149,A[1:150,],main="Game of Life, 500x500\nInitial Density Test with dead borders",xlab="time",ylab="density")
fields:::image.plot(y=rnd_threshold,x=0:149,B[1:150,],main="Game of Life, 500x500\nInitial Density Test with border wrapping",ylab="density",xlab="time")
text(20,.5,"Initial densities from .01 and .96 all\nconverge to somewhere between .38-.39\n.97+ crash to 0\n.6 to .96 crash too, but recover",pos=4)


# does edge-awareness explain this?
PerimeterRatio <- function(side){
	(side*4-4)/(side^2)
}
# http://s19.postimage.org/9dp2lqqb7/Percent_Edge_Aware.png

plot(0:249,sapply(seq(500,2,-2),PerimeterRatio),type='l',
		main="proportion edge-aware by time on 500x500 GOL board",
		xlab="time",ylab="proportion",ylim=c(0,1))