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