Multifractals in ecology using R - Day 4

Posgraduate course, University of Maringá, 2013

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

Leonardo A. Saravia

Percolation

Exercise 1

Infection

Exercise 2