Replicating NetLogo Fire model
[This article was first published on R snippets, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
While preparing for the new semester I have started reimplementing standard NetLogo examples in R. The first is Fire model.Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
The simulation in R is presented here:
# Forest matrix trees encoding:
# 3 – green, 2 – burning, 1 – burnt, 0 – no tree
simulate <- function (forest.density = 0.6, size = 251) {
forest <- matrix(3 * rbinom(size ^ 2, 1, forest.density),
nrow = size)
forest[1,] <- 2
while (sum(forest == 2) > 0) {
list.red <- which(forest == 2)
# find spots touching burning trees
ignite.x <- 1 + floor((list.red – 1) / size)
ignite.x <- rep(ignite.x, each = 4) +
rep(c(1, –1, 0,0), length(ignite.x))
ignite.y <- 1 + ((list.red – 1) %% size)
ignite.y <- rep(ignite.y, each = 4) +
rep(c(0, 0, 1,–1), length(ignite.y))
include <- (ignite.x > 0) & (ignite.x <= size) &
(ignite.y > 0) & (ignite.y <= size)
ignite.n <- (ignite.x[include] – 1) * size +
ignite.y[include]
# ignite green trees
is.green <- (forest[ignite.n] == 3)
forest[ignite.n[is.green]] <- 2
# trees burn for one period
forest[list.red] <- 1
}
return(1 – sum(forest[-1,] == 3) / sum(forest[-1,] > 0))
}
The only thing left was to check that the results produced by NetLogo and R implementations are identical. I have simulated both models 1000 times for forest density equal to 60%. Here is the plot comparing percent burned density estimates and Kolmogorov-Smirnov test p-value.
It can be concluded that for the parametrization I have selected models produce similar results.
To leave a comment for the author, please follow the link and comment on their blog: R snippets.
R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.