STAT 360 Lab 19 & 20: Spatial Simulations (Parts 2 and 3)
Name: Rebecca Lewis
Set Up Work Space
Load R Libraries
library(rmarkdown)
library(knitr)Set the Seed
set.seed(33)The Bak et al. (1990) Fire Model
Background:
The following lattice model was originally proposed by Bak et al. (1990; Phys Lett A 147:297-300). Consider an 'n' by 'n' grid of cells where each cell is either empty (value of 0), occupied by a healthy tree (value of 1), or occupied by a burning tree (value of 2). At each time step, (1) empty cells are recolonized by healthy trees at a rate of 'p', (2) burning trees spread fire to healthy trees in neighboring cells with a probability of 'r', and (3) burning trees burn down and leave behind empty cells.Constructing the fireSim Function
Exercise 1:
Use function() to construct a function and assign it to the object fireSim. This function should have the following attributes and perform the following computations:
(1) Take in as arguments the objects n, p, r, and t, and display with default
values of 20, 0.5, 0.85, 10, and TRUE, respectively
(2) Initialize a blank data frame with t rows and three columns with the names
"Empty", "Healthy", and "Burning" and assign it to the object totalTrees
(3) Initialize two blank matrices with n rows and n columns each and assign
them to the objects lattice.now and lattice.next
(4) Fill all cells in lattice.now with healthy trees
(5) Use sample() to transform a single cell in lattice.now from a healthy tree
to a burning tree
(6) Use for() to construct a loop that will repeat the following computations t
times:
(a) Set all values in lattice.next to NA
(b) Use for() to construct a nested for loop that will cycle through
every cell in lattice.now and perform the following computations for
each cell:
(i) Use if() to determine if the cell in lattice.now is empty; if
it is, use if() and runif() to determine if it should be
recolonized; if is should, assign a healthy tree to that cell
in lattice.next, else assign an empty cell there
(ii) Use if() to determine if the cell in lattice.now is occupied
by a healthy tree; if it is, assign a healthy tree to that
cell in lattice.next
(c) Use for() to construct a nested for loop that will cycle through
every cell in lattice.now and perform the following computations for
each cell:
(i) Use if() to determine if the cell in lattice.now is occupied
by a burning tree; if it is, perform the following:
(A) Use if() to determine if the cell to the right of it is
occupied by a healthy tree; if it is, use if() and
runif() to determine if that healthy tree should catch
on fire; if it should, assign a burning tree to that
cell in lattice.next
(B) Use if() to determine if the cell to below it is
occupied by a healthy tree; if is is,use if() and
runif() to determine if that healthy tree should
catch on fire; if it should, assign a burning tree to
that cell in lattice.next
(C) Use if() to determine if the cell to the left of it is
occupied by a healthy tree; if it is, use if() and
runif() to determine if that healthy tree should catch
on fire; if it should, assign a burning tree to that
cell in lattice.next
(D) Use if() to determine if the cell above it is occupied
by a healthy tree; if it is, use if() and runif() to
determine if that healthy tree should catch on fire; if
it should, assign a burning tree to that cell in
lattice.next
(ii) Use if() to determine if the cell in lattice.now is occupied
by a burning tree; if it is, assign an empty cell to that
cell in lattice.next
(d) Use if() to determine if the argument display is equal to TRUE; if it
is, use image(t(lattice.now), col = c("gray", "green", "red)) to
display lattice.now and Sys.sleep(0.25) to pause the loop before
continuing to its next iteration
(e) Set lattice.now equal to lattice.next
(f) Use sum() to calculate the total number of 0s, 1s, and 2s in
lattice.now and assign them to the corresponding row in totalTrees
(7) Use return() to return totalTrees once the function is complete
Hint 1: You can use the segments of code designed during Lab 18 for many of
these computations
Hint 2: The segments of code designed during Lab 18 were specific to the 1st
row and 1st column of lattice.now - you will need to update those
references to correspond to the row and column of the nested for loop
Hint 3: Fire cannot spread beyond the boundaries of the lattice - the function
if() will come in handy to determine if not all 4 nearest neighbors
should be checked when determining transitions to burning trees
Hint 4: It will help to run segments of the code individually (which should not
be knitted) before combining everything into a single functionfireSim<-function(p=.5,r=.85,t=10,display=TRUE, n=20)
{
TotalTrees<-data.frame(Empty=rep(NA,t), Healthy=rep(NA, t),Burning=rep(NA,t))
lattice.now<-matrix(nrow = n, ncol = n)
lattice.next<-matrix(nrow = n, ncol = n)
lattice.now[1:n,1:n]<-1
row<-sample(1:nrow(lattice.now),1,replace=FALSE)
col<-sample(1:ncol(lattice.now),1,replace=FALSE)
lattice.now[row,col]<-2
for(u in 1:t)
{
lattice.next[1:n,1:n]<-NA
for(i in 1:n)
{
for(j in 1:n)
{
if(lattice.now[i,j]==0)
{
if(runif(n = 1, min = 0, max = 1) < p)
{
lattice.next[i,j]=1
}
else{lattice.next[i,j]=0}
}
if(lattice.now[i,j]==1)
{
lattice.next[i,j]=1
}
}
}
for(i in 1:n)
{
for(j in 1:n)
{
if(lattice.now[i,j]==2 && j!=20)
{
if(lattice.now[i,(j+1)]==1)
{
if(runif(n = 1, min = 0, max = 1) < r)
{
lattice.next[i,(j+1)]=2
}
}
}
if(lattice.now[i,j]==2 && i!=20)
{
if(lattice.now[(i+1),j]==1)
{
if(runif(n = 1, min = 0, max = 1) < r)
{
lattice.next[(i+1),j]=2
}
}
}
if(lattice.now[i,j]==2 && j!=1)
{
if(lattice.now[i,(j-1)]==1)
{
if(runif(n = 1, min = 0, max = 1) < r)
{
lattice.next[i,(j-1)]=2
}
}
}
if(lattice.now[i,j]==2&& i!=1)
{
if(lattice.now[(i-1),j]==1)
{
if(runif(n = 1, min = 0, max = 1) < r)
{
lattice.next[(i-1),j]=2
}
}
}
if(lattice.now[i,j]==2)
{
lattice.next[i,j]=0
}
}
}
if(display==TRUE)
{
image(t(lattice.now),col=c("gray","green","red"))
Sys.sleep(0.25)
}
lattice.now<-lattice.next
TotalTrees$Empty[u]<-sum(lattice.now==0)
TotalTrees$Healthy[u]<-sum(lattice.now==1)
TotalTrees$Burning[u]<-sum(lattice.now==2)
}
return(TotalTrees)
}Using the fireSim Function
Exercise 2:
Use fireSim() to project the spread of a fire on a 20 by 20 grid with p = 0.5 and r = 0.8 for t = 100 time steps and without displaying the lattice and assign it to the object fire1.
fire1<-fireSim(p<-.5,r<-.8,t<-100,n<-20, display=FALSE)Exercise 3:
Use plot(), lines(), and fire1 to generate a time plot of the number of empty cells and healthy and burning trees over the t time steps of the simulation. Hint: the colors should match the color scheme set up in fireSim.
fire1<-fireSim(p<-.5,r<-.8,t<-100,n<-20, display=FALSE)
plot(fire1$Empty, ylim=c(0,400), main="Forest Fire", xlab="Time Steps", ylab="Number of Trees")
lines(fire1$Empty, type="l", col="gray")
lines(fire1$Healthy, type="l", col="green")
lines(fire1$Burning, type="l", col="red")
Exercise 4:
Estimate how long it takes the system to settle into a stable (but oscillating) pattern.
Around 20 times. Exercise 5:
Calculate the maximum number of burning trees observed during this simulation and display the result.
max(fire1$Burning)## [1] 81Exercise 6:
Use replicate() and fireSim() to simulate the spread of fire using the same parameters as above 100 times, calculate the average maximum number of burning trees observed during the 100 simulations, and display the result.
BurnBaby<-mean(replicate(100, max(fireSim(p<-.5,r<-.8,t<-100,n<-20, display=FALSE)$Burning)))
BurnBaby## [1] 82.99Exercise 7:
Use replicate() and fireSim() to simulate the spread of fire using the same parameters as above 100 times, calculate the average total number of burning trees observed during the 100 simulations, and display the result.
Burn<-mean(replicate(100, sum(fireSim(p<-.5,r<-.8,t<-100,n<-20, display=FALSE)$Burning)))
Burn## [1] 5952.07Exercise 8:
Use for() to construct a loop that will cycle through every value of r between 0.5 and 1.0 (in intervals of 0.05 and while holding all else constant), use replicate and fireSim() to record the average total number of burning trees observed during 100 simulations to the object totalBurned, and display the result. Hint: all other simulation parameters should remain unchanged.
totalBurned <- c(.5,.55,.6,.65,.7,.75,.8,.85,.9,.95,1)
for(i in seq(.5, 1, by=.05))
{
if(i==.5)
{
totalBurned[1]<-mean(replicate(100,sum(fireSim(r=i,display=FALSE,p=.05,n=20,t=10)$Burning)))
}
if(i==.55)
{
totalBurned[2]<-mean(replicate(100,sum(fireSim(r=i,display=FALSE,p=.05,n=20,t=10)$Burning)))
}
if(i==0.60)
{
totalBurned[3]<-mean(replicate(100,sum(fireSim(r=i,display=FALSE,p=.05,n=20,t=10)$Burning)))
}
if(i==.65)
{
totalBurned[4]<-mean(replicate(100,sum(fireSim(r=i,display=FALSE,p=.05,n=20,t=10)$Burning)))
}
if(i==.7)
{
totalBurned[5]<-mean(replicate(100,sum(fireSim(r=i,display=FALSE,p=.05,n=20,t=10)$Burning)))
}
if(i==.75)
{
totalBurned[6]<-mean(replicate(100,sum(fireSim(r=i,display=FALSE,p=.05,n=20,t=10)$Burning)))
}
if(i==.8)
{
totalBurned[7]<-mean(replicate(100,sum(fireSim(r=i,display=FALSE,p=.05,n=20,t=10)$Burning)))
}
totalBurned[8]<-mean(replicate(100,sum(fireSim(r=.85,display=FALSE,p=.05,n=20,t=10)$Burning)))
if(i==.9)
{
totalBurned[9]<-mean(replicate(100,sum(fireSim(r=i,display=FALSE,p=.05,n=20,t=10)$Burning)))
}
if(i==.95)
{
totalBurned[10]<-mean(replicate(100,sum(fireSim(r=i,display=FALSE,p=.05,n=20,t=10)$Burning)))
}
if(i==1)
{
totalBurned[11]<-mean(replicate(100,sum(fireSim(r=i,display=FALSE,p=.05,n=20,t=10)$Burning)))
}
}
totalBurned## [1] 30.09 50.70 68.82 89.21 102.38 117.26 129.07 135.11 143.02 146.89
## [11] 150.70Exercise 9:
Use plot() to generate a scatterplot of the average total number of trees burned as a function of r and use it to determine which value of r between 0.5 and 1.0 (in intervals of 0.05) results, on average, in the smallest total number of trees burned. Hint: this is an example of a brute force approach to finding an optimum set of parameters.
plot(totalBurned, main="Burning Trees", xlab=" The value of r", ylab="Number of Burnt Trees",xaxt = 'n')
axis(1, at=1:11, c(.5,.55,.6,.65,.7,.75,.8,.85,.9,.95,1))
The smallest value of burning trees on average is the r value .5Adjusting the fireSim Function
Exercise 10:
Use function() to re-create fireSim such that the fire can now spread to the 8 nearest neighbors instead of the 4 nearest neighbors and assign it to the object fireSim2.
fireSim2<-function(p=.5,r=.85,t=10,display=TRUE, n=20)
{
TotalTrees<-data.frame(Empty=rep(NA,t), Healthy=rep(NA, t),Burning=rep(NA,t))
lattice.now<-matrix(nrow = n, ncol = n)
lattice.next<-matrix(nrow = n, ncol = n)
lattice.now[1:n,1:n]<-1
row<-sample(1:nrow(lattice.now),1,replace=FALSE)
col<-sample(1:ncol(lattice.now),1,replace=FALSE)
lattice.now[row,col]<-2
for(u in 1:t)
{
lattice.next[1:n,1:n]<-NA
for(i in 1:n)
{
for(j in 1:n)
{
if(lattice.now[i,j]==0)
{
if(runif(n = 1, min = 0, max = 1) < p)
{
lattice.next[i,j]=1
}
else{lattice.next[i,j]=0}
}
if(lattice.now[i,j]==1)
{
lattice.next[i,j]=1
}
}
}
for(i in 1:n)
{
for(j in 1:n)
{
if(lattice.now[i,j]==2 && j!=20)
{
if(lattice.now[i,(j+1)]==1)
{
if(runif(n = 1, min = 0, max = 1) < r)
{
lattice.next[i,(j+1)]=2
}
}
}
if(lattice.now[i,j]==2 && i!=20)
{
if(lattice.now[(i+1),j]==1)
{
if(runif(n = 1, min = 0, max = 1) < r)
{
lattice.next[(i+1),j]=2
}
}
}
if(lattice.now[i,j]==2 && i!=20 && j!=20)
{
if(lattice.now[(i+1),(j+1)]==1)
{
if(runif(n = 1, min = 0, max = 1) < r)
{
lattice.next[(i+1),(j+1)]=2
}
}
}
if(lattice.now[i,j]==2 && i!=1 && j!=20)
{
if(lattice.now[(i-1),(j+1)]==1)
{
if(runif(n = 1, min = 0, max = 1) < r)
{
lattice.next[(i-1),(j+1)]=2
}
}
}
if(lattice.now[i,j]==2 && i!=20 && j!=1)
{
if(lattice.now[(i+1),(j-1)]==1)
{
if(runif(n = 1, min = 0, max = 1) < r)
{
lattice.next[(i+1),(j-1)]=2
}
}
}
if(lattice.now[i,j]==2 && i!=1 && j!=1)
{
if(lattice.now[(i-1),(j-1)]==1)
{
if(runif(n = 1, min = 0, max = 1) < r)
{
lattice.next[(i-1),(j-1)]=2
}
}
}
if(lattice.now[i,j]==2 && j!=1)
{
if(lattice.now[i,(j-1)]==1)
{
if(runif(n = 1, min = 0, max = 1) < r)
{
lattice.next[i,(j-1)]=2
}
}
}
if(lattice.now[i,j]==2&& i!=1)
{
if(lattice.now[(i-1),j]==1)
{
if(runif(n = 1, min = 0, max = 1) < r)
{
lattice.next[(i-1),j]=2
}
}
}
if(lattice.now[i,j]==2)
{
lattice.next[i,j]=0
}
}
}
if(display==TRUE)
{
image(t(lattice.now),col=c("gray","green","red"))
Sys.sleep(0.25)
}
lattice.now<-lattice.next
TotalTrees$Empty[u]<-sum(lattice.now==0)
TotalTrees$Healthy[u]<-sum(lattice.now==1)
TotalTrees$Burning[u]<-sum(lattice.now==2)
}
return(TotalTrees)
}Exercise 11
Use replicate() and fireSim2() to simulate the spread of fire using the same parameters (as those used in Exercise 7) 100 times, calculate the average maximum number of burning trees observed during the 100 simulations, and display the result.
BurnBabyBurn<-mean(replicate(100, max(fireSim2(p<-.5,r<-.8,t<-100,n<-20, display=FALSE)$Burning)))
BurnBabyBurn## [1] 113.56