【问题标题】:How can I shorten the runtime of for loops and if statements in R while using igraph for forest fire simulations如何在使用 igraph 进行森林火灾模拟时缩短 R 中 for 循环和 if 语句的运行时间
【发布时间】:2020-03-15 21:03:30
【问题描述】:

我在 R 中模拟森林火灾,必须使用 igraph 包。我的代码目前可以工作,但速度极慢。我通读了矢量化 for 循环或使用 seq_along 或将条件放在循环之外的方法。我无法弄清楚如何在我的特定代码中使用这些解决方案。至于我的代码的描述:我正在模拟森林火灾,其中我循环了 21 个不同的百分比,代表一个空白顶点变成一棵树的可能性(0 到 1,间隔为 0.05)。在每个循环中,我都在运行 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")

【问题讨论】:

  • 我喜欢这种锻炼方式。你确实有很多嵌套循环,并且在你内心深处,你在每个子集的基础上运行neighbours()。您可能可以省略循环并让 igraph 直接在邻居组上工作。首先运行一些基准测试,检查你的代码在你的机器上哪里慢,这样你就可以把精力集中在正确的地方。从事情最慢的地方开始。下周晚些时候,当我有空闲时间时,我将纯粹出于灵感来解决这个问题。我非常乐观地认为它可以显着加速。祝你好运。
  • 你确定V(H)$color[V(G)[color == "black"][i]] &lt;- "green"?你真的想把黑树变成绿色吗?

标签: r loops if-statement simulation igraph


【解决方案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(sum( V(G)$color=="red" ), 0, 1)

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

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

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

当您查找消耗大量时间的迭代时,您会发现您的脚本会因检查燃烧树木的邻居和燃烧的邻居而受到很大影响。

引入新行为以促进优化:

我的优化解决方案是为已经蔓延的火灾引入一种新颜色:“橙色”。由于所有具有燃烧邻居的树木在每个时间步长内都会着火,因此模拟不需要检查在前一个时间步长之前咳嗽的树木的邻居。这显着减少了 adjacent_vertices() 执行的邻居测试的数量,该函数将在 p=.05 上运行 20*100*50*270 次左右。那是一百万个邻居检查!如果我们不需要检查已经点亮所有邻居的黄色树的邻居,我们可以节省大量 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 ")

【讨论】:

    猜你喜欢
    • 1970-01-01
    • 2016-02-28
    • 2021-06-19
    • 1970-01-01
    • 2021-10-12
    • 1970-01-01
    • 2017-06-02
    • 2011-07-10
    • 1970-01-01
    相关资源
    最近更新 更多