Model Thinking- Aggregation- Game of Life- new cycle function with wrapping instead of dead borders

GOLwrapping.r

### FIXED!!

# It turns out that the assumptions you make for border have a huge impact on the evolution of a 
# Game of Life board. Our first implementations (Vadim and myself both) were to assume dead cells on the 
# borders. Here I've included a new function where edge cells have neigbors on the opposite side of the board
# i.e. the board wraps around both top to bottom and right to left, which is a shape called a taurus.

# this function runs about twice as fast as the previosu best-performer from Vadim and gives VERY different 
# results (just looking at one output, measurement, the proportion of cells alive at any given iteration).

# main finding: when borders wrap around, the system stabilizes faster and maintains a higher density
# when borders are dead (0s), the system asymptotically stabilizes, since roving cells formations will
# often get absorbed into the edges (die off). 

# given this observation, the two kinds of systems shouldbe studied SEPARATELY, and it's hard to speak of
# GOL dynamics in general without first identifying what the border assumptions are.

# first, to compare, here are Vadim's functions again:

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)
}

# and now here is my GOL function where borders wrap around (just 1 function)
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] <- 1
	return(M)
}

# simulate to see if it makes a difference:
set.seed(1)
M1 <- matrix(sample(c(1,0),size=10000,replace=TRUE,prob=c(.2,.8)),ncol=100)
prop_on100 <- vector(length=1001)
prop_on100[1] <- sum(M1)/10000
for (i in 1:1000){
	M1 <- GOLborderwrap(M1)
	prop_on100[i+1] <- sum(M1)/10000
}

set.seed(1)
M2 <- matrix(sample(c(1,0),size=10000,replace=TRUE,prob=c(.2,.8)),ncol=100)
prop_on100_2 <- vector(length=1001)
prop_on100_2[1] <- sum(M2)/10000
for (i in 1:1000){
	M2 <- life_cycle(M2)
	prop_on100_2[i+1] <- sum(M2)/10000
}

set.seed(1)
M1 <- matrix(sample(c(1,0),size=1000000,replace=TRUE,prob=c(.2,.8)),ncol=1000)
prop_on1000 <- vector(length=1001)
prop_on1000[1] <- sum(M1)/1000000
for (i in 1:1000){
	M1 <- GOLborderwrap(M1)
	prop_on1000[i+1] <- sum(M1)/1000000
}

set.seed(1)
M2 <- matrix(sample(c(1,0),size=1000000,replace=TRUE,prob=c(.2,.8)),ncol=1000)
prop_on1000_2 <- vector(length=1001)
prop_on1000_2[1] <- sum(M2)/1000000
for (i in 1:1000){
	M2 <- life_cycle(M2)
	prop_on1000_2[i+1] <- sum(M2)/1000000
}

plot(1:1001,prop_on100,type='l',col="blue",ylim=c(0,1),main="It makes a really big difference\nwhat border assumptions you make!")
lines(1:1001,prop_on100_2,col="red")
lines(1:1001,prop_on1000,col="blue",lty=2)
lines(1:1001,prop_on1000_2,col="red",lty=2)
legend("topleft",lty=c(1,1,2,2),col=c("blue","red","blue","red"),legend=c("wrap100","dead100","wrap1000","dead1000"))



# that's my two cents. So what is preferred, dead or wrapping borders?!