fire.R

st <- proc.time()

# given const
N        = 50     # grid size
PercSteps= 20     # number of percolation values
pstep    = 1/(PercSteps) # density step
reps     = 50     # number of repetitions
Perc     = matrix(0, PercSteps+1, 3)

# allMap <- vector('list', PercSteps+1)
# allMapN <- 1

# APAN: r 是 1:PercSteps+1 還是 1:(PercSteps+1) ? 我改成後者 1:(PercSteps+1)
for (r in 1:(PercSteps+1)){

  # cat("\tr = ", r)

  p      = pstep*(r-1)
  Time   = vector(mode = "integer", reps) # 用 integer 即可

  for (k in 1:reps){ # k 為完整燒到結束, 共完整燒 reps 次

    full     = matrix(0L, N+2, N+2)  # define running matrix (to set the boundary of the grid with 0, so that the fire wont go out )
    # .Random.seed <- 1:1000
    full[ 2:(N+1) , 2:(N+1) ] <- rbinom(N**2, 1, p) # 中間填二元亂數
    treeLoc <- which(full == 1, arr.in=T) # 樹的位置, 後面有用
    full[2, (1:N)+1] [full[2, (1:N)+1] == 1] <- 2L # 第一排有火
    # allMap[[allMapN]] <- full
    count = 0L;         # count number of iterations

    # while (sum(full==2)>0){
    while (any(full == 2)){

      count = count + 1L
      y <- which(full==2L, arr.in=TRUE)   # find trees on fire
      size = nrow(y)

      #### 以下這部份才是重點，也就是算火怎麼傳染的過程
      #### 你原本的演算法好懂直覺, 但太多 if() 以及有些不需要的轉換
      #### 這也是會很慢的原因
      #### 我的方法概念上大概是「每個時刻都找到所有會燒的樹，把它加上1」
      #### 經過比較，運算結果和你的版本完全相同，但可快至少一倍

      # 往上座標 y - matrix( c(1,0), ncol=2, nrow=size , byrow=T )
      # 往下座標 y + matrix( c(1,0), ncol=2, nrow=size , byrow=T )
      # 往左座標 y - matrix( c(0,1), ncol=2, nrow=size , byrow=T )
      # 往右座標 y + matrix( c(0,1), ncol=2, nrow=size , byrow=T )
      fireTreeLoc <- rbind(
        y,
        y - matrix( c(1,0), ncol=2, nrow=size , byrow=T ),
        y + matrix( c(1,0), ncol=2, nrow=size , byrow=T ),
        y - matrix( c(0,1), ncol=2, nrow=size , byrow=T ),
        y + matrix( c(0,1), ncol=2, nrow=size , byrow=T )
      ) # 目前和下次會燒起來的座標, 但仍包括沒樹的位置
      fireTreeLoc.abs <- (fireTreeLoc[,2] - 1)*(N+2) + fireTreeLoc[,1]
      treeLoc.abs <- (treeLoc[,2] - 1)*(N+2) + treeLoc[,1]
      addOneLoc <- treeLoc.abs[treeLoc.abs %in% fireTreeLoc.abs]
      # addOneLoc <- fireTreeLoc.abs[fireTreeLoc.abs %in% treeLoc.abs] #同義

      fireLayer <- matrix(0L, (N+2), (N+2)) ; fireLayer[addOneLoc] <- 1L
      full <- full + fireLayer
      # full[full > 3] <- 3 # 不必要, 只是把 > 3 的值硬改成 3, 但對運算結果沒作用

      # allMap[[allMapN]] <- full
      # allMapN <- allMapN + 1
    }
    Time[k] = count
  }

  Mean_Time = mean(Time)
  Error_Time = sd(Time)
  Perc[r,1] = p
  Perc[r,2] = Mean_Time
  Perc[r,3] = Error_Time
}

plot(Perc[,1],Perc[,2]/N)
print(Perc)
print(proc.time() - st)




# setwd("/tmp")
# png(file="fire%05d.png")
# for(i in 10000:11000){
#   map <- data.frame(z=c(allMap[[i]]), y=rep(1:52,52) , x=as.numeric(gl(52,52)))
#   print(ggplot(map, aes(x=x,y=y,fill=z)) + geom_tile(aes(fill=z)) + scale_fill_gradient(low="green", high="red"))
# }
# dev.off()