mercredi 20 novembre 2019

How can I shorten the runtime of for loops and if statements in R while using igraph for forest fire simulations

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