Model Thinking- Aggregation- Cellular Automata- Experiments in R

CellularAutomata.R

# the code executes cellular automata 2^(2^3) experiments under 2^3 neighbor rules. Where each rule determines whether a cell turns on or off (TRUE-FALSE). Code is semi-efficient in turns of vectorizing processes, but there is plenty of room for improvement. 

# first the rule matrix 8 rows (rules) by 3 columns (positions relative to middle)
ruleMatrix <- matrix(0,ncol=3,nrow=8)
ruleMatrix[2,3] <- 1
ruleMatrix[3,2] <- 1
ruleMatrix[4,2:3] <- 1
ruleMatrix[5,1] <- 1
ruleMatrix[6,c(1,3)] <- 1
ruleMatrix[7,1:2] <- 1
ruleMatrix[8,1:3] <- 1

# the reference matrix for neighbor rules.
# each of these 8 rows is selected according to
# the rows of 'BigRuleMatrix', defined later
image(x=-1:1,y=1:8,t(ruleMatrix),axes=FALSE,ylab="",xlab="rule",asp=1,main="rule matrix")
axis(1,at=-1:1,lwd=0)
text(-1.5,1:8,labels=8:1,xpd=TRUE,pos=2)

# functions:

# determines the 3 cell values (1 left, same, and 1 right)
# given a row (vector) of values and an index number
linear3Neighbors <- function(ind,vec){
	inds <- c(ind-1,ind,ind+1)
	inds[inds==0] <- NA
	vec[inds]
}

# given the 3 cell values (from 'linear3Neighbors()'), does the middle cell live or die? (on-off)
cellOn <- function(ind,onRules,vec){
	n3 <- linear3Neighbors(ind,vec)
	any(apply(onRules,1,function(x,n3) {all(n3 == x) },n3 = n3),na.rm=TRUE)
}

# a function to make the 256 binary rule switches:
BinComb <- function(nrules=8){
	BigRuleMatrix <- matrix(nrow=2^nrules,ncol=nrules)
	powers <- 2^(0:(nrules-1))
	for (i in 1:nrules){
		BigRuleMatrix[, i] <- rep(c(rep(1,powers[(8-i+1)]),rep(0,powers[(8-i+1)])),powers[i])
	}
	return(BigRuleMatrix)
}

# all of the binary rule combinations
BigRuleMatrix <- BinComb()
image(BinComb()) # take a look: better image in final pdf

# define a container to hold results BigResultsArray
BigResultsArray <- array(dim=c(201,220,256))
print(object.size(BigResultsArray),units="Mb") # 42 Mb, not bad

# this sets the random number seed so that the same results
# can be reproduced: the only ranomness enters with the 
# generation of the first row. Since each rule set shares the
# same starting conditions, everything else is determined
# by the first row. Different starting conditions will generate
# different results. Trying this exhaustive experiment out under
# a large set of starting conditions would be a strong motivation
# to parallelize this code, since a single run takes so long..

# each row is 220 elements. where the first row is FALSE for the first and last 100 entries, and
# randomly generated for the middle 20 entries (more likely FALSE than TRUE)

set.seed(1)
N <- 200 # the number of iterations to do...
# first row will be identical for all runs: this may take 10-30 minutes,
# depending on your machine... The z loop can be parallelized as well,
# but I didn't take that step.
row1 <-  c(rep(FALSE,100),sample(c(TRUE,FALSE),20,prob=c(.1,.9),replace=TRUE),rep(FALSE,100))
for (z in 1:256){
	onRules <- matrix(as.logical(ruleMatrix[as.logical(BigRuleMatrix[z,]),]),ncol=3)
	Results <- matrix(ncol=220,nrow=(N+1))
	Results[1,] <-  row1
	
	# the i loop 'must' be done with a loop- cannot be vectorized...
	for (i in 1:N){
		Results[i+1,] <- sapply(1:ncol(Results),cellOn,onRules=onRules,vec=Results[i,])
	}
	# save results in z level of array
	BigResultsArray[,,z] <- Results
	gc() # clears garbage memory, not necessary
}

# BigResultsArray contains everything we need to produce graphical output:

pdf("WolframWithRandomStarts.pdf")
# page 1: the 8x3 rule matrix:
image(x=-1:1,y=1:8,t(ruleMatrix),axes=FALSE,ylab="",xlab="rule",asp=1,main="rule matrix")
axis(1,at=-1:1,lwd=0)
text(-1.5,1:8,labels=8:1,xpd=TRUE,pos=2)
text(-5.8,4,"This is the rule matrix used\nthroughout, where 1-8 are all\npotential rules. Light yellow = 'on'\nand red = 'off'; -1 is\none position to the left and 1\nis one to the right\nRules 1-8 are switched on and\noff according to the matrix on the\nnext page",pos=4,xpd=TRUE)
# page 2: the rule matrix switches:
image(x=1:8,y=1:256,t(BinComb())[,256:1],main="Yellow = 'on', Red = 'off'",asp=.25,axes=FALSE,xlab="",ylab="")
for (i in 1:8){
	segments(9,seq(256-i+1,1,by=-8)-.7,10.3+i,seq(256-i+1,1,by=-8)-.7,col=paste(gray(.5),50,sep=""),lwd=.5)
	text(10+i,seq(256-i+1,1,by=-8)-.5,labels=seq(i,256,by=8),pos=4,cex=.34,xpd=TRUE)
}
text(-40,125,"Each row of this matrix is the\n'switch' pattern for the rule matrix\nshown on the previous page\nthere are 256 possible switch patterns\n(2^8)",pos=4,xpd=TRUE)
axis(1,at=1:8,labels=NA)
text(1:8,-7,1:8,xpd=TRUE,cex=.8)
text(-40,60,"In following, each image starts with the \nsame randomly generated row and\niterates downward(200 times)\naccording to rules from the previous page\nactivated by the switches shown here and\nnumbered accordingly",pos=4,xpd=TRUE)
# page 3-257
# each results page. This exercise can be redone with new random starts. But
# it only makes senese if each rule set shares the same starting conditions. Raster keeps the
# resulting file a bit smaller.
for (i in 1:256){
image(t(BigResultsArray[,,i])[,201:1],useRaster=TRUE,main=paste("rule configuration =",i),axes=FALSE)
}
dev.off()
cat("the big pdf is located here:",getwd())