Model Thinking- Tipping Points- Percolation

Percolation.R

# Author: triffe
###############################################################################

# this code will try to implement an example percolation model, where density and matrix size are parameters

# the main function also produces neat statistics along the way for post-analysis
# dead borders. either rook or queen neighborhoods work.

# um, name means percolation, with dead borders
PercIteration <- function(M, neighborhood.type = "rook"){
	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(0, rbind(0, M, 0), 0)
	
	# how many neighbors?:
	# left side
	M.out <- M
	if (neighborhood.type == "queen"){
		M.out[M.aug[1:nr, 1:nc] > 0]  		    <- 1
		M.out[M.aug[2:(nr+1), 1:nc] > 0]  	    <- 1
		M.out[M.aug[3:(nr+2), 1:nc] > 0]  	    <- 1
			
					# middle:  		
		M.out[M.aug[1:nr,2:(nc+1)] > 0]  	    <- 1
		M.out[M.aug[3:(nr+2),2:(nc+1)] > 0]  	    <- 1
			
					# right:
		M.out[M.aug[1:nr,3:(nc+2)] > 0]  	    <- 1
		M.out[M.aug[2:(nr+1),3:(nc+2)] > 0]  	    <- 1
		M.out[M.aug[3:(nr+2),3:(nc+2)] > 0]  	    <- 1
	}
	if (neighborhood.type == "rook"){
		M.out[M.aug[2:(nr+1), 1:nc] > 0] 	    <- 1
		M.out[M.aug[1:nr,2:(nc+1)] > 0] 	    <- 1
		M.out[M.aug[3:(nr+2),2:(nc+1)] > 0]	    <- 1
		M.out[M.aug[2:(nr+1),3:(nc+2)] > 0]	    <- 1
	}
	
	# birth and death rules:
	M[M.out > 0 & !is.na(M)] <- 1
	
	return(M)
}

PercMatrix <- function(ncol = 100, nrow = 100, dens = .59){
	M <- matrix(ncol = ncol, nrow = nrow)
	M[sample(1:(ncol * nrow), size = floor(ncol * nrow * dens), replace = FALSE)] <- 0
	# let's say the fire starts in the first row:
	M[1, M[1, ] == 0] <- 1
	M
}

Percolate <- function(M, image = TRUE, sleep = 0, col = c("green3","firebrick"), PercIteration = PercIteration){
	if (image){
		image(M, col = col, useRaster = TRUE, main = "t = 0")
	}
	# starts
	i 		<- 0
	nb2 	        <- sum(M, na.rm = TRUE)
	N 		<- length(M[!is.na(M)])
	s.dens 	        <- N / length(M)
	any.new         <- TRUE
	
	# stats to take along the way, preallocate overly large matrix 
	# and chop off NAs at end
	stats <- matrix(nrow = .3*prod(dim(M)), ncol = 4)
	colnames(stats) <- c("it","p.inf","n.inf.it","trav")
	# stats for starting matrix
	stats[1, 1] <- i # iteration
	stats[1, 2] <- sum(M, na.rm = TRUE) / N
	stats[1, 3] <- nb2 
	stats[1, 4] <- any(M[ nrow(M), ] == 1)
	# it = iteration
	# p.inf = percent infiltrated:
	# n.inf.it = number cells infiltrated this iteration
	# trav = has the matrix been traversed yet? (will indicate time at which traversed, although inf can continue)
	
	# while loop! watch out! hehe
	while (any.new){
		i <- i + 1
		nb1 <- nb2
		M <- PercIteration(M)
		if (image){
			Sys.sleep(sleep)
			image(M, col = col, useRaster = TRUE, main = paste("t =", i))
		}
		nb2 <- sum(M, na.rm = TRUE)
		any.new <- nb2 > nb1
		
		stats[i+1, 1] <- i # iteration
		stats[i+1, 2] <- sum(M, na.rm = TRUE) / N
		stats[i+1, 3] <- nb2 - nb1
		stats[i+1, 4] <- any(M[, ncol(M)] == 1)
	}
	stats <- stats[!is.na(stats[, 1]), ]
	out <- list(stats = stats, t.traverse = ifelse(stats[nrow(stats), 4], min(stats[stats[, 4], 1]), NA), N = N, s.dens = s.dens)
}

# example of how it works:
stats <- Percolate(PercMatrix(ncol = 100, nrow = 100, dens = .59), image = TRUE)
# stats is a list with 4 components:
names(stats)
# take a look at prop infiltrated over time:
plot(stats$stats[, 1], stats$stats[, 2], type = 'l', xlab = "time", ylab = "proportion infiltrated",
		main = paste("Percolation Model\nRook neighborhood\nstarting density of", stats$s.dens))

# ----------------------------------------------------------------------------------------------
# just for an idea of the kinds of things you can do with the output data:

# infiltration hazard over time:
x <- stats$stats[, 1]
# N is the number of cells that could be infiltrated (non NAs):
N <- stats$stats$N
# S is the proportion remaining non-infiltrated
lx <- 1 - stats$stats[, 2]
# dx is the probability density function of being infiltrated- it need not have a regular shape:
dx <- c(1-lx[1], -diff(lx))
plot(x, dx, type = 'l')
# these are exposures:
Lx <- lx + .5 * dx
# mx is the hazard rate:
mx <- dx / Lx
plot(x, mx, type = 'l')
# average time spent non-iniltrated of the cells that will eventually be infiltrated:
(e0inf <- sum(Lx - Lx[length(Lx)]))

# ----------------------------------------------------------------------------------------------
# And if you wanted to reproduce the curve shown in lectures, where the density tipping point is made obvious:
# this is a super-fine excercise, so do with caution- i.e. if you have a big awesome computer no prob
# but if you're on a laptop, then reduce the number of densities to test over and reduce the number of trials
# per density: (This test running on a remote server at the time of this writing, will post results soon)

#densities <- seq(0,1,length.out=1001)
## we'll do each density 100 times!
#reps <- 1000
#library(parallel)
#cl <- makeCluster(detectCores())
#results <- parSapply(cl, densities, function(x, reps, Percolate, PercMatrix, PercIteration){
#			 replicate(reps, !is.na(Percolate(PercMatrix(ncol = 100, nrow = 100, dens = x), image = FALSE, PercIteration = PercIteration)$t.traverse))
#		}, reps = reps, Percolate = Percolate, PercMatrix = PercMatrix, PercIteration = PercIteration)
#stopCluster(cl)
#colnames(results) <- densities
#rownames(results) <- 1:reps
#save(results, file = "PercolationTipping.Point.Rdata")
#
# p.traverse <- colSums(results) / reps
# plot(densities, p.traverse, type ='l', xlab = "density", ylab = "probability full traverse", main = # "Probility of making full traverse by density\naverage of 1000 reps per density in intervals of .001", )

# to do the same on a single core (lots of time, don't be so ambitious):

densities <- seq(0,1,length.out=101)
# we'll do each density 100 times!
reps <- 200

results <- sapply(densities, function(x, reps, Percolate, PercMatrix, PercIteration){
			replicate(reps, !is.na(Percolate(PercMatrix(ncol = 100, nrow = 100, dens = x), image = FALSE, PercIteration = PercIteration)$t.traverse))
		}, reps = reps, Percolate = Percolate, PercMatrix = PercMatrix, PercIteration = PercIteration)
colnames(results) <- densities
rownames(results) <- 1:reps
save(results, file = "PercolationTipping.Point.Rdata")