I am simulating forest fires in R and have to use the igraph package. My code currently works but is extremely slow. I read through ways of vectorizing my for loops or using seq_along or putting conditionals outside my loops. I was unable to figure out how to use these solutions in my specific code. As for the description of my code: I am simulating forest fires where I loop through 21 different percentages representing the likelihood of a blank vertex becoming a tree (0 through 1 by .05 intervals). In each of these loops I am running 100 full forest fires. Each forest fire is comprised of 50 time steps. In each time step, I check which vertices of my igraph need to be changed to empty, tree, and fire. For the specific problem I am working on, I am tracking the largest number of trees on fire during each forest fire so that I can later generate a graph of the average maximum fire for the 21 different percentages. Any tips on how to speed up this code would be much appreciated.
OG <- graph.lattice(c(30,30))
V(OG)$color <- "black"
total.burning.tree.max <- matrix(nrow = 21, ncol = 100)
for (p in seq(0, 1, .05)) {
for (x in 1:100) {
fire.start <- sample(900, 1)
tree.start <- sample(900, (900*.7))
G <- OG
V(G)$color[tree.start] <- "green"
V(G)$color[fire.start] <- "red"
current.burning.tree.max <- 1
H <- G
for (h in 1:50) {
if (length(V(G)[color == "red"]) > current.burning.tree.max) {
current.burning.tree.max <- length(V(G)[color == "red"])
}
for (i in 1:length(V(G)[color == "black"])) {
if (runif(1) <= p) {
V(H)$color[V(G)[color == "black"][i]] <- "green"
}
}
if (length(V(G)[color == "red"]) > 0) {
for (d in 1:length(V(G)[color == "red"])) {
V(H)$color[V(G)[color == "red"][d]] <- "black"
potential.fires <- neighbors(G, V(G)[color == "red"][d])
for (z in 1:length(potential.fires)) {
if (V(G)$color[potential.fires[z]] == "green") {
V(H)$color[potential.fires[z]] <- "red"
}
}
}
}
G <- H
}
total.burning.tree.max[(p*20), x] <- current.burning.tree.max
print(current.burning.tree.max)
}
}
burn.numbers <- c()
for (c in 1:21) {
burn.numbers[c] <- average(total.burning.tree.max[c, ])
}
plot(burn.graph, type = "l")
Aucun commentaire:
Enregistrer un commentaire