优化代码的一般注意事项:
首先,您的代码充满了嵌套循环,其中每个模拟循环遍历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 ")