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

Comments