在使用igraph进行森林火灾模拟时，如何缩短R中for循环和if语句的运行时间

Question

我在R中模拟森林火灾，必须使用igraph包。我的代码目前可以工作，但速度非常慢。我通读了向量化for循环的方法，或者使用seq_along的方法，或者将条件句放在循环之外的方法。我不知道如何在我的特定代码中使用这些解决方案。至于我的代码的描述:我正在模拟森林火灾，我循环遍历21个不同的百分比，表示空白顶点变成树的可能性(通过.05间隔从0到1)。在每个循环中，我运行了100个完整的森林火灾。每一次森林火灾由50个时间步长组成。在每个时间步中，我检查我的igraph的哪些顶点需要更改为空、树和火。对于我正在处理的特定问题，我正在跟踪每次森林火灾期间着火的最大树木数量，以便稍后可以生成21个不同百分比的平均最大火灾的图表。任何关于如何加速这段代码的技巧都将不胜感激。

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")

Answer 1

关于优化代码的一般说明：

首先，您的代码充满了嵌套循环，每个模拟循环遍历igraph中的节点以更改值。这不是一个好主意，因为igraph更快。

例如，像这样在给定颜色的所有节点上循环：

for (i in 1:length(V(G)[color == "red"])) {
  V(H)$color[V(G)[color == "red"][i]] <- "black"
}

最好是存储节点的子集，并使用它一次性进行所有更改：

V(G)[ V(G)$color=="red" ] <- "black"

还要注意，您不需要将runif(1, p)放在循环中，但是如果让runif()输出一个向量，则可以执行任意数量的概率比较:runif($color==(V(G)runif(1, p)“red”)，0，1) <= p

当您不需要变量或igraph节点属性的实际值时，可以考虑汇总布尔值：

sum(V(G)$color=="red") == length( V(G)$color[ V(G)$color  =="red" ] )

在您的示例中，通常在运行常规模拟或在igraph中运行模拟时，计算速度取决于模拟中的动态。例如，我下面的脚本执行的时间步长要快得多，几乎没有着火的树。当函数adjacent_vertices()被指示返回mode="total"时，它在这里是一个明显的时间盗贼。然而，这个函数应该比你自己循环更快。

当你寻找耗费大量时间的迭代时，你会发现你的脚本在检查燃烧的邻居和燃烧的邻居时会受到很大的影响。

引入新的行为以促进优化：

我的优化解决方案是引入一种新的颜色：“橙色”，用于已经蔓延的火灾。由于在每个时间步中所有邻居着火的树都着火了，因此模拟不需要检查在前一个时间步之前着火的树的邻居。这大大减少了在p=.05上运行20*100*50*270次左右的adjacent_vertices()的邻居测试次数。这是一百万的邻居检查就在那里！如果我们不需要检查已经拥有所有邻居的黄树的邻居，我们就节省了大量的CPU周期。

​

​

我希望我已经提供了一些很好的通用提示。在上面的脚本旁边，我希望下面的脚本可以用于教学目的。

在下面的脚本中，我更改了存储模拟数据的方式，以及模拟中我可能无法理解的一个函数。下面的p现在说明了每个时间步扑灭燃烧的树的概率，而燃烧树的邻居肯定会在下一个时间步着火(就像它们在模拟中一样)。

每个级别的p都绘制了一个示例图。

还请注意，可以通过删除runif()来稍微优化设置新树着火的行，该行允许您更改相邻树着火的单独概率的值。

tree_fires <-  potential_fires[  runif(length(potential_fires), 0, 1) <= FIRE_PROBABILITY  ]

一如既往的优化。把你的努力花在他们有价值的地方！与迁移到橙色树来简化adjacent_vertices()的工作相比，删除tree_fires的runif()可能只节省了大约百万分之一的时间。

关于你的方法的注意事项：

我也做过类似的社交网络中死亡传播的模拟。你把最初的火放在哪里很重要。在一次迭代中着火的树的最大数量被你的森林的墙覆盖了很多。这将导致在p假设的每个级别内您的测量结果有更大的变化。我非常建议您转移到一个模型，该模型将最初的火灾放置在您的森林中间。我已经为此包含了配置变量。

建议摘要：

library("igraph")
# Configurations
PROB_LEVELS <- 20            # How many probability levels?
FOREEST_SIMULATIONS <- 100   # How many simulations shouls occur for each probability level?
TIMESTEPS <- 50              # How many iterations shouls fires spread for in each simulation?
FIRE_PROBABILITY <- 1        # How likely is it that an adjacent tree will catch fire? (Lower values decrease speed of fire spreading)
FIXED_STARTING_POINT <- TRUE # Should the fire begin at the same place always?
PLAYGROUND <- 30             # The size of the forest (higher values decrease likelyhood of hiting foret-walls)
FOREST_DENSITY <- .7         # The percentage of nodes that are trees in an unburnt forest. (higher values facilitates spread of fire)

# 900 trees
OG <- graph.lattice(c(PLAYGROUND, PLAYGROUND))
V(OG)$color <- "gray"
# Store simulation results in a list instead.
stat <- lapply(1:PROB_LEVELS, function(x) rep(NA,FOREEST_SIMULATIONS))

plotforest <- function(graph){plot(graph, vertex.label=NA, vertex.size=5, layout=layout_on_grid(graph) )}

# Make dimulations using these probabilities
for (p in 1:PROB_LEVELS/PROB_LEVELS) {
  cat("p =",p)

  for (x in 1:FOREEST_SIMULATIONS) {
    # Each iteration have different random configurations of forests with a fixed tree-density
    G <- OG
    V(G)$color[ sample(PLAYGROUND^2, (PLAYGROUND^2 * FOREST_DENSITY )) ] <- "green"
    # Firees could start at random tree or in the "middle"
    if(FIXED_STARTING_POINT){
      V(G)$color[ round(PLAYGROUND^2/2)-(PLAYGROUND/2) ] <- "red" }
    else{
      V(G)$color[ sample(PLAYGROUND^2, 1) ] <- "red" }


    # Collect simulation data over time-steps during which the fire spreads
    burning_tree_max <- 1
    for(h in 1:TIMESTEPS){
      # Put out trees that are on fire using probability `p`
      # This replaces your loop for (i in 1:length(V(G)[color == "red"])) {}
      trees_on_fire <- V(G)[ V(G)$color=="red" ] # make this subset only once per iteration. Store it. You could use %in% c('red','orange' )
      if(length(trees_on_fire) == 0){break;print(h)} # Abort time-steps if there are no more contageous fires.
      V(G)$color[ trees_on_fire[ runif(length(trees_on_fire), 0, 1) <= p ] ] <- "black"

      # Set neighboring trees of burning trees on fire (only green trees can catch fire)
      # This replaces your nested loop staring with for (d in 1:length(V(G)[color == "red"])) { }
      last_egnited <-  V(G)$color=="red"
      potential_fires <- adjacent_vertices(G, last_egnited, mode="total")
      potential_fires <- unique(unlist(potential_fires))
      potential_fires
      tree_fires <-  potential_fires[  runif(length(potential_fires), 0, 1) <= FIRE_PROBABILITY  ]
      # Store last time-step's burning trees as orange, and egnite new neighbors.
      V(G)$color[last_egnited] <- "orange"
      V(G)$color[tree_fires][V(G)$color[tree_fires] == "green"] <- "red" # Set all green subsetted neighbors of flaming treas on fire at once
      # No orange tree can have a green neighbour!

      # Track maximum number of trees on fire.
      burning_tree_max <- max(burning_tree_max, sum(V(G)$color=="red") )
    }

    # store simulation results as sum of currently burning trees
    stat[[p*PROB_LEVELS]][x] <- burning_tree_max

  }
  cat(": averaging", round(mean(stat[[p*PROB_LEVELS]], na.rm=T),1), "trees.", fill=T)
  plotforest(G)

}


# Plot the simulation results
plot(sapply(stat, function(x) mean(x)), type="l",
     ylab="Maximum number of trees on fire", xlab=NA,
     main="Snapshot of fires during a simulation",
     sub="50 time-cycles ona 30x30 sized forest ")