# Percolation in R - Day 4

## Percolation

• First we will generate a random forest: set each position of a matrix to 1 or 0 with some probability $p$

• First we need to fill a matrix with numbers

		rm(list=ls())

dim_rf <- 10

matrix(0,dim_rf,dim_rf)

matrix(1:dim_rf,dim_rf)

matrix(1:(dim_rf*dim_rf),dim_rf)

• what we really need are random numbers
		runif(10)

rf <- matrix(runif(dim_rf*dim_rf),ncol=dim_rf)

• now we need to decide which will become 0 or 1 with some probability, the best way to do that is using a function
		set_tree <- function(elem,prob) {
if(elem>prob)
elem <- 0
else
elem <- 1
}

set_tree(rf,0.5)

• the last line don’t works because we need to apply the condition to each site individually
		rf <-apply(rf,1:2,set_tree,0.5)

• So we set a matrix of 1’s and 0’s with probability 0.5, we can plot it
		image(rf,asp=1,axes=F,col=c("white","black"))

• Now we can put all together and make a function that does all things
		generate_rf <- function(sideDim,p)
{
rf <- matrix(runif(sideDim*sideDim),ncol=sideDim)
rf <-apply(rf,1:2,set_tree,p)
image(rf,asp=1,axes=F,col=c("white","black"))
}

• We can see how the matrix is filled when we change $p$
		generate_rf(100,0.1)

generate_rf(100,0.5)

generate_rf(100,0.59)


• The next thing is to burn the trees, we forget to return the matrix with the random forest
		generate_rf <- function(sideDim,p)
{
rf <- matrix(runif(sideDim*sideDim),ncol=sideDim)
rf <-apply(rf,1:2,set_tree,p)
image(rf,asp=1,axes=F,col=c("white","black"))
return(rf)
}

rf <- generate_rf(10,0.51)

• Now we need a function to explore the neighborhood and fire the actual site if one on the adjacent sites is fired. We can try to use a loop.
		rf[1,] <- 2

for(j in 1:ncol(rf))
{

if( rf[1,j]==2 & rf[2,j]==1)
rf[2,j] <- 2
}

• But we have to do that for all the rows of the matrix
		for(i in 2:nrow(rf))
for(j in 1:ncol(rf))
{
if( rf[i-1,j]==2 & rf[i,j]==1)
rf[i,j] <- 2
}

image(rf,asp=1,axes=F,col=c("white","black","red"))

• something is wrong, we need to test all the neighborhood
		for(i in 2:nrow(rf))
for(j in 2:(ncol(rf)-1))
{
if( (rf[i-1,j]==2 || rf[i-1,j-1]==2 || rf[i-1,j+1]==2) && rf[i,j]==1)
rf[i,j] <- 2
}

• So now we are ready to build the function
		fire_rf <- function(rf) {
dimi=nrow(rf)
dimj=ncol(rf)-1
rf[1,] <- 2
for(i in 2:dimi)
for(j in 2:dimj)
{
if((rf[i-1,j]==2 || rf[i-1,j-1]==2 || rf[i-1,j+1]==2) && (rf[i,j]==1))
rf[i,j] <- 2
}
return(rf)
}

rf <- generate_rf(100,0.1)

rf1 <- fire_rf(rf)

image(rf1,asp=1,axes=F,col=c("white","black","red"))

• Now to finish we have to count the number of burned sites, a simple function will do, but first we test the commands
   		bur <-0

for(i in 2:nrow(rf1))
for(j in 1:ncol(rf1))
{
if(rf1[i,j]==2 )
bur <- bur + 1
}

bur

• next we create the function
		countBurned <- function(rf)
{
bur <-0
dimi=nrow(rf)
dimj=ncol(rf)

for(i in 2:dimi)
for(j in 1:dimj)
{
if(rf1[i,j]==2 )
bur <- bur + 1
}
return(bur/((dimi-1)*dimj))
}


## Exercise 1

• Make a plot of the proportion of burned sites versus the probability $p$ of the trees

## Infection

• We can do this in a similar way, first we need a matrix but now one dimension will be the time
		dim_in <- 10
time_in <- 20

inf<-matrix(0,time_in,dim_in)

• At time 1 we need to infect some sites to have a start
		inf[1,] <- ifelse(runif(dim_in)>0.5,1,0)

• Now we have to propagate the infection, there are two possibilities: contagion with probability $\lambda$ or recuperation with probability $\mu$
		lambda = 0.5

mu= 0.5

for(i in 1:(time_in-1))
for(j in 2:(dim_in-1))
{
if(inf[i,j]==0){

if(inf[i,j-1]==1){
inf[i+1,j] <- ifelse(runif(1)<=lambda,1,0)
}
else if(inf[i,j+1]==1 ){
inf[i+1,j] <- ifelse(runif(1)<=lambda,1,0)
}
}
else
{

inf[i+1,j] <- ifelse(runif(1)<=mu,0,1)
}
}

image(inf,asp=1,axes=F,col=c("grey","brown"))

• Then we put all in a function
		simulate_inf <- function(lambda,mu,dim_in,time_in){

inf<-matrix(0,time_in,dim_in)
inf[1,] <- ifelse(runif(dim_in)>0.5,1,0)
for(i in 1:(time_in-1))
for(j in 2:(dim_in-1))
{
if(inf[i,j]==0){

if(inf[i,j-1]==1){
inf[i+1,j] <- ifelse(runif(1)<=lambda,1,0)
}
else if(inf[i,j+1]==1 ){
inf[i+1,j] <- ifelse(runif(1)<=lambda,1,0)
}
}
else
{
inf[i+1,j] <- ifelse(runif(1)<=mu,0,1)
}
}
image(inf,asp=1,axes=F,col=c("grey","brown"))
}

simulate_inf(.4,.4,50,100)

• We can try with different parameters and see what happens at $\lambda > \mu$ or $\lambda < \mu$

## Exercise 2

• Build the plot with a fixed $\mu$ of the probability of propagation versus $\lambda$

• Estimate the fractal dimension and the multifractal spectrum of the infection

Previous