Skip to content

Instantly share code, notes, and snippets.

@timriffe

timriffe/PercolationSim.R

Created March 21, 2012 10:26
Show Gist options
  • Save timriffe/2146034 to your computer and use it in GitHub Desktop.
Save timriffe/2146034 to your computer and use it in GitHub Desktop.
Download ZIP
Model Thinking- Tipping Points- Percolation- Simulation to find tipping point
Raw
PercolationSim.R
# Author: Tim Riffe
# this code (almost entirely commented out because it would take so long to re-run [however, you have the
# option!]) shows how the below results were derived. This percolation code uses a stripped-down version of
# the previous percolation code here: https://gist.github.com/2044583 - it only produces 1 result: does the
# infiltration reach the other side of the board or not, using the rook definition. The code was reduced to
# be as concise as possible and then compiled, to run as fast as possible (though I'm sure there's an even
# faster way)
# the flow: we want to try out the same percolation model for each starting density from 0 to 1 in increments
# of .001 on a 100x100 board. Each density we repeat 300 times. Since the output is TRUE/FALSE (traversed or
# didn't) we can take the mean of the length 300 vector, and this gives us an estimate of the probability
# for that density. If you skip to the bottom, I pasted the results vector using dput(), and the only active
# code here is the plot of the results.
densities <- seq(0,1,length.out=1001)
# we'll do each density 100 times!
reps <- 300
# the simulation is commented out because it seriously takes an entire day or more
# to run!
#library(parallel)
#cl <- makeCluster(detectCores())
#results <- parSapply(cl, seq(0,1,length.out=1001), function(x, reps){
# PercIteration <- function(M, nc, nr){
# # 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
# 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){
# nc <- ncol(M)
# nr <- nrow(M)
# # starts
# i <- 0
# nb2 <- sum(M, na.rm = TRUE)
# any.new <- TRUE
# # while loop! watch out! hehe
# while (any.new){
# i <- i + 1
# nb1 <- nb2
# M <- PercIteration(M, nc = nc, nr = nr)
# nb2 <- sum(M, na.rm = TRUE)
# any.new <- nb2 > nb1
# }
# return(any(M[nrow(M),] == 1, na.rm = TRUE))
# }
#
# PercIteration <- compiler:::cmpfun(PercIteration)
# Percolate <- compiler:::cmpfun(Percolate)
# mean(replicate(reps, Percolate(PercMatrix(ncol = 100, nrow = 100, dens = x))))
# }, reps = 300)
#stopCluster(cl)
#
#colnames(results) <- densities
#rownames(results) <- 1:reps
#save(results, file = "PercolationTipping.Point.Rdata")
##
#p.traverse <- colSums(results) / reps
#file.choose()
#load("C:\\Users\\Administrador\\Desktop\\PercolationTipping.Point.Rdata")
# results from dput(p.traverse)
p.traverse <- c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.00333333333333333,
0.00333333333333333, 0.01, 0, 0.00333333333333333, 0, 0.00333333333333333,
0.00333333333333333, 0.00666666666666667, 0.01, 0.01, 0.00666666666666667,
0.00666666666666667, 0.02, 0.01, 0.01, 0.01, 0.0166666666666667,
0.02, 0.0366666666666667, 0.0266666666666667, 0.0466666666666667,
0.04, 0.0333333333333333, 0.0466666666666667, 0.0666666666666667,
0.0266666666666667, 0.0733333333333333, 0.0733333333333333, 0.1,
0.06, 0.123333333333333, 0.113333333333333, 0.133333333333333,
0.166666666666667, 0.16, 0.163333333333333, 0.18, 0.196666666666667,
0.23, 0.206666666666667, 0.25, 0.263333333333333, 0.303333333333333,
0.323333333333333, 0.38, 0.383333333333333, 0.393333333333333,
0.426666666666667, 0.45, 0.43, 0.473333333333333, 0.523333333333333,
0.503333333333333, 0.536666666666667, 0.603333333333333, 0.57,
0.646666666666667, 0.636666666666667, 0.7, 0.713333333333333,
0.703333333333333, 0.743333333333333, 0.816666666666667, 0.776666666666667,
0.783333333333333, 0.803333333333333, 0.87, 0.86, 0.89, 0.866666666666667,
0.883333333333333, 0.89, 0.906666666666667, 0.923333333333333,
0.923333333333333, 0.953333333333333, 0.966666666666667, 0.953333333333333,
0.963333333333333, 0.936666666666667, 0.963333333333333, 0.966666666666667,
0.993333333333333, 0.97, 0.986666666666667, 0.983333333333333,
0.973333333333333, 0.993333333333333, 0.983333333333333, 0.996666666666667,
0.986666666666667, 1, 0.996666666666667, 0.993333333333333, 0.996666666666667,
0.993333333333333, 1, 1, 0.996666666666667, 1, 1, 1, 1, 0.996666666666667,
1, 0.996666666666667, 0.996666666666667, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1)
densities <- seq(0,1,length.out=1001)
plot(densities, p.traverse, type ='l', xlab = "density", ylab = "probability full traverse", main = "Probility of making full traverse by density\naverage of 300 reps per density in intervals of .001", )
# this is the density closest to a .5 probability of traversing
densities[which.min(abs(p.traverse-.5))] # 0.594
# a couple of tippiness functions, re: the last posted lecture from Scotte- these are not generalized to
# multiple outcomes: these will only work for 2 outcomes, where x is a vector of probabilties, and the other
# possible outcome is just 1-x
# diversity:
divI <- function(x){
(x ^ 2 + (1 - x) ^ 2) ^ -1
}
# entropy
entI <- function(x){
-x * log(x, base = 2)
}
# what's the diversity index function?
plot(seq(0,1,by=.001), divI(p.traverse), type='l', main = "diversity index of percolation model", xlab = "starting density", ylab = "diversity")
# what's the entropy index function?
plot(seq(0,1,by=.001), entI(p.traverse), type='l', main = "entropy index of percolation model", xlab = "starting density", ylab = "entropy")
# End--
# suggestion: do same with queen neighborhood: ought to yield similar curve but shifted left. Then
# where would the tipping point be?
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment