Percolation Algorithm In R

I’m casually following along with a few of the lessons in Princeton’s Algorithms, Part 1 on Coursera. The course itself is taught in Java, but I thought completing some of the lessons in R would improve my R skills. The first activity I’ve completed is a percolation simulation using weighted quick-union with path compression. I’m not actually sure if this was one of the assignments, but it seemed cool so I thought I’d give it a try.

Percolation is meant to model physical or social systems where some type of blockage occurs. The algorithm taught in the class is to begin with full blockage and to randomly unblock one element until there is a path from one side of the system to another. The simplest method of doing this is to imagine a blocked grid. The initial blocked system is shown below.

percolation start

My algorithm removes one blockage per half second and this is modeled by turning the black box white. A totally percolated system is shown below. Notice there is now a path through the blockage – by starting at the bottom, following the open (white) squares up and to the left and then diagonally up and to the right.

percolation end

At the moment percolation is completed a message is printed to the console.

percolation output

The code to initialize the board and percolated system and to run the clearing algorithm is shown below.

####################################################################
# James McCammon
# 21 Feb 2015
# Percolation Algorithm Functions Version 1
####################################################################
# Load libraries
library(shape)
source("perc_functions.R")
 
# Define variables
nrows = 10
ncols = 10
nsites = nrows*ncols
closedSites = 1:nsites
numneighbors = 4
display = NULL
#display = matrix(data=0, nrow=nrows, ncol=ncols)
VTN = nsites + 1
VBN = nsites + 2
 
# Initailize board
board = initBoard()
initDisplay()
initVTN()
initVBN()
 
# Run main
main()

Created by Pretty R at inside-R.org

Most of the work is done in the set of functions below. This can definitely be improved, for example by incorporating an arbitrary list of neighbors into each site object rather than assuming the system is a square grid and determining neighbors by looking in the four quadrants around each site. I could also think harder about memory usage. Right now I just threw everything into the global environment and didn’t think too hard about memory usage or efficiency. Using the tracemem() function from the pryr shows there are a few times the board object is copied. I’ll probably play around with things and continue to improve them.

####################################################################
# James McCammon
# 21 Feb 2015
# Percolation Algorithm Main Version 1
####################################################################
# Define main function
main = function() {
  sleeptime = 0.5
 
  while(!connected(VTN, VBN)) {
    randSite = getRandSite()
    board[[randSite]]$open = TRUE
    #location = trans2display(randSite)
    #write2display(location)
    filledrectangle(mid=c(display[[randSite]]$xmid, 
                          display[[randSite]]$ymid), 
                    display[[randSite]]$wx, 
                    display[[randSite]]$wy, 
                    col='white', 
                    lwd=1, 
                    lcol='white')
 
    neighbors = getNeighbors(randSite)
 
    for(i in 1:numneighbors) {
      if(is.null(neighbors[[i]])) next
      if(!isOpen(neighbors[[i]])) next
      join(randSite, neighbors[[i]])
    } 
    Sys.sleep(sleeptime)
  }
  print("Percolation completed")
  p = 1 - (length(closedSites)/nsites)  
  return()
}
 
 
#Initalize the display
initDisplay = function() {
  display <<- as.vector(1:nsites, mode='list')
  windows(5, 5)
  plot.new()
  xmid = 0
  ymid = 1
  wx = .1
  wy = .1
  index = 1
 
  for(row in 1:nrows) {
    for(col in 1:ncols) {
      filledrectangle(mid=c(xmid, ymid), wx=wx, wy=wy, col='black', lwd=1, lcol='white')
      display[[index]] <<- list('xmid'=xmid, 'ymid'=ymid, 'wx'=wx, 'wy'=wy)
      index = index + 1
      xmid = xmid + .1
    }
    xmid = 0
    ymid = ymid - .1
  }  
}
 
 
# Gets root 
getRoot = function(node) {
  while(node != board[[node]]$root) {
    board[[node]]$root <<- board[[board[[node]]$root]]$root
    node = board[[node]]$root
  }
  return(node)
}
 
 
# Checks if two notes are connected
connected = function(node1, node2) {
  return(getRoot(node1) == getRoot(node2))
}
 
 
# Initalize board
initBoard = function() {
  board = as.vector(1:nsites, mode='list')
  board = lapply(board, function(i) list("index"=i, "size"=1, "root"=i, "open"=FALSE))
  return(board)
}
 
 
# Initalize virtural top node
initVTN = function() {
  board[[VTN]] <<- list("index"=VTN, "size"=(ncols + 1), "root"=VTN, "open"=FALSE)
  lastrow = nsites - ncols + 1
  for(i in 1:ncols)  {
    board[[i]]$root <<- board[[VTN]]$root
  }  
}
 
 
# Initalize virtual bottom node
initVBN = function() {
  board[[VBN]] <<- list("index"=VBN, "size"=(ncols + 1), "root"=VBN, "open"=FALSE)
  lastrow = nsites - ncols + 1
  for(i in lastrow:nsites) {
    board[[i]]$root <<- board[[VBN]]$root
  }  
}
 
 
# Choose random site
getRandSite = function() {
  randSite = sample(closedSites, size=1)
  closedSites <<- closedSites[-which(closedSites==randSite)]
  return(randSite)
}
 
 
# Checks if site is open
isOpen = function(site) {
  return(!(site%in%closedSites))
}
 
 
# Get neighbors on a square board
getNeighbors = function(index) {
  top = index - ncols
  if(top < 1) top = NULL
 
  bottom = index + ncols
  if(bottom > nsites) bottom = NULL
 
  left = index - 1
  if(!(left%%ncols)) left = NULL
 
  right = index + 1
  if(right%%ncols == 1) right = NULL
 
  return(list("tn"=top, "rn"=right, "bn"=bottom, "ln"=left))
}
 
 
# Join two nodes
join = function(node1, node2) {
 
  r1 = getRoot(node1)
  r2 = getRoot(node2)
  s1 = board[[r1]]$size
  s2 = board[[r2]]$size
 
  if(s1 > s2) {
    board[[r2]]$root <<- board[[r1]]$root
    board[[r1]]$size <<- s1 + s2
  }
 
  else {
    board[[r1]]$root <<- board[[r2]]$root
    board[[r2]]$size <<- s2 + s1 
  } 
}
 
 
####################################################################
# Functions for visualizing output with matrix
####################################################################
# Translate location of site to matrix
trans2display = function(index) {
  col = index%%ncols
  if(!col) col=ncols
  row = ceiling(index/nrows)
  return(list("row"=row, "col"=col))  
}
 
 
# Update the matrix display
write2display = function(transObject) {
  display[transObject$row, transObject$col] <<- 1  
}

Created by Pretty R at inside-R.org

Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s