# Multifractals in Ecology Using 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
1 2 3 4 5 6 7 8 9

`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
1 2 3

`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
1 2 3 4 5 6 7 8

`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
`1`

`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
`1`

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

- Now we can put all together and make a function that does all things
1 2 3 4 5 6

`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\)
1 2 3 4 5

`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
1 2 3 4 5 6 7 8 9

`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.
1 2 3 4 5 6 7 8

`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
1 2 3 4 5 6 7 8 9

`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
1 2 3 4 5 6

`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
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18

`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
1 2 3 4 5 6 7 8 9 10

`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
1 2 3 4 5 6 7 8 9 10 11 12 13 14

`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
1 2 3 4

`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
`1`

`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\)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24

`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
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25

`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