Data Analysis Perspectives, Computing in R

flee circus banner

What The Flea Circus Can Teach Us About Computing In R

By Richard Herrington, consultant, UNT University Information Technology

Recently, fellow R&SS consultant Jon Starkweather brought a computing problem to my attention – originally this was  posted on the Project Euler project site – Problem 213.

The problem description is as follows: Flea Circus (Problem 213):  A 30×30 grid of squares contains 900 fleas, initially one flea per square. When a bell is rung, each flea jumps to an adjacent square at random (usually four possibilities, except for fleas on the edge of the grid or at the corners).   What is the expected number of unoccupied squares after 50 rings of the bell?  Give your answer rounded to six decimal places. 

Jon's approach to the following problem can be found at:  http://bayes.acs.unt.edu:8083/BayesContent/class/Jon/R_SC/Module12/Fleas_001.R.  At first glance,  many R users might approach the problem by using looping constructs (e.g. for,  while,  repeat).   However, R programmers also know that the key to efficient programming in R in using vectorized operations,  rather than using for loops to iterate over the indices of vector and matrix objects.  We explore the following:  Can vectorized approaches be used in the Flea Circus problem?  Does the structure of the problem allow other efficient programming constructs to be used? This brief note explores using the following approaches in R:    i) recursive function calls;  and  ii) parallel implementation of the for-loop across multiple CPUs;  and iii) byte-compiling of user created functions; - in programming the Flea Circus problem in R.  According to Jon, his R code, takes approximately 45 minutes for about 250 replications of the 50-step iteration of 900 fleas jumping on a 30 by 30 grid (using for and while loops,  with no multiple CPU parallelization, or byte compiling of functions calls). I will explore how much speed-up we can achieve in a Monte Carlo implementation of the Flea Circus problem, by using these techniques in conjunction. We will find that we can perform 250 replications in approximately six minutes on a single CPU (with NO parallelization and using recursive functions that are  byte-compiled). Parallelization across 8 cores takes less than a minute for 250 iterations. This is approximately 1.6 seconds for all 50-steps of iteration for all 900 fleas.  Parallelization of 2000 iterations (250 iterations  * 8 core) takes approximately 6.4 minutes on 8 core -  roughly the same amount of time with 8 times  the number of iterations. Published solutions to the  Project Euler problems can be found here. The published solution with 6 digits of accuracy) is 330.721154.

Overview of the R Scripts and the General Approach Taken

The approach taken here is to keep track of the locations of the particles (fleas) as they move around on the 30 by 30 grid of cell locations - move.incr.  The Flea Circus example is an example of a discrete markov chain process - a 30 by 30 matrix of probabilities that are iterated 50 times to produce a final matrix of probabilities,  that is finally used to calculate the expected values of the empty cells at the end of 50 steps of iteration.  This type of iteration process (when the output from the function becomes the input to the same function) suggests the use of a function that calls itself until a criterion is met, to end the iterations – a recursive function.  

We will create a recursively called function - flea.jump - such that at the end of a sequence of 50 steps of recursive calls to itself,  the resulting expected value of the empty cells is returned to the function labled flea.circus.repetition.  The function flea.circus.repitition  is recursively called until the number of pre-specified replications has occured.   Some of the run-time speed decreases of the script are a result of using an apply function to the move list using the recursively called function produce.valid.move (note: using a for-loop here would slow down the function down considerably, hence the use of the apply function). The produce.valid.move function is called recursively untill a valid move is produced for a single current flea location.   This updated move is returned to the apply function within the flea.jump function. The final output of that apply function is an updated move list for ALL 900 fleas.  Within the function flea.jump,  a  table of probabilities is generated for the current locations of the fleas - after jumping. The complement of these probabilities on step i,  are P_i = (1 - p) and are normalized  to sum to 1, and then they are multiplied by the corresponding matrix of probabilities from the last iteration (i.e. P_i-1). Lastly, this resulting multiiplication of the two probabilites (P_i-1  *  P_i ) is again normalized to sum to 1, which is used to produce the current matrix of probabilities for empty cells. From this 30 by 30 matrix of probabilities,  the sum of the probabilities that correspond to empty cells on the current step,  are multiplied by the constant 900 cells,  to produce the expected value of the empty.cell count on the current iteration. The empy.cell count and flea.mat.current.prob probability values  are then saved into the flea.mat.list structure to be used within the next iteration of flea.jump. The function flea.circus.repetition is recursively called to replicate the entire 50 steps of the 900 fleas,  a number of pre-specified times - number.of.repetitions.  The expected empty.cell count is accumulated over iterations, and the replicated average is returned once the flea.circus.repetition  is finished.  

Unfortunately, R is limited to no more than a couple thousand nested recursions,  therefore,  we need to implement an outer controlling function that uses a loop to replicate the recursively called replications of the function flea.circus.repetition.  This function - monte.carlo.rep.func can be used with R packages which are designed to utilize multiple CPUs on a single machine, or across multiple machines with multiple CPUs (e.g. R packages:  foreach, parallel, doMC). This parallization of the for-loops across multiple CPUs will allow us to gain further time reductions, for a  larger numbers of iterations of the function monte.carlo.rep.func.

A Note About The Precision of Recursive Arithmetic Computations

The original problem statement requires that the results should be reported with six digits of precision. Iterated arithmetic computations are susceptable to accumulating round-off error,  at each step of the calcuation. Inacurate results will be produced, after a certain number of iterations, if the accuracy of the computer's internal numeric representations are not sufficient.  Additionally, Monte-Carlo variation across replications will produce varying estimates of the final calculations.   If the accuracy of the computer's floating-point-unit is not very accurate,  and Monte-Carlo variation is large,  then it will be difficult (maybe impossible) to achieve six digits of accuracy. While we don't explore these issues here, it might be beneficial to use R packages that are designed to produce arbitrarily accurate precision for arithmetic calculations in R (e.g. R package Rmpfr).   

High Performance Computing at UNT

If the use of parallelization in your research project(s) sounds appealing to you, and you don't have access to high performance mutlicore machines (more than four CPUs), then you might take a look at University Information Technology's high performance computing cluster. As stated on their website:  "The cluster systems are intended for computationally intensive LINUX-capable software. Systems may be available for use by any  faculty or currently enrolled qualified student at UNT, with additional requirements to use some resources. Student access must be sponsored by a faculty member."  The UIT's HPC cluster includes R (and R packages) as part of UIT's HPC maintained suite of software

R Function Declarations

Function  produce.valid.move

# Load package to byte-compile functions
require(compiler)

produce.valid.move<-function(location)
  {location.new<-location
   axis.to.move<-sample(c("x", "y"), 1)
   step.incr<-sample(c(1, -1), 1)
   if(axis.to.move=="x") location.new[1]<-location.new[1] + step.incr
   if(axis.to.move=="y") location.new[2]<-location.new[2] + step.incr
   if(location.new[1]>=1 & location.new[1]<=30 & location.new[2]>=1 & location.new[2]<=30) 
      return(location.new)
   produce.valid.move(location)}

# Byte-compile the function
produce.valid.move<-cmpfun(produce.valid.move)

Function flea.jump

flea.jump<-function(flea.mat.list)
    {current.step<-flea.mat.list$current.step+1
     empty.cells.old<-empty.cells
     flea.mat.current.prob.old<-flea.mat.list$flea.mat
     move.incr<-t(apply(flea.mat.list$flea.mat.location, 1, produce.valid.move))
     flea.mat.current.freq<-matrix(table(move.incr[,1], move.incr[,2]), nrow=30, ncol=30)
     flea.mat.current.prob.temp<-flea.mat.current.freq + 1
     flea.mat.current.prob.temp<-flea.mat.current.prob.temp/sum(flea.mat.current.prob.temp)
     flea.mat.current.prob<- (flea.mat.current.prob.old) * (1-flea.mat.current.prob.temp)
     flea.mat.current.prob<-flea.mat.current.prob/sum(flea.mat.current.prob)
     empty.cells<-sum(flea.mat.current.prob[flea.mat.current.freq==0])*900
     flea.mat<-flea.mat.current.prob
     flea.mat.list<-list(empty.cells=empty.cells,flea.mat=flea.mat.current.prob,
                         flea.mat.location=move.incr, current.step=current.step)
     if(current.step==50)
       {
       return(flea.mat.list)
       }
     flea.jump(flea.mat.list)}

# Byte-compile the function
flea.jump<-cmpfun(flea.jump)

Function flea.circus.repetition

flea.circus.repetition<-function(flea.mat.list, number.of.repetitions, current.iteration, 
    empty.cells.total)
   {gc()
    cat(current.iteration,"\n")
    empty.cells.old <-empty.cells.total
    flea.mat<-matrix(1, nrow=30, ncol=30)
    flea.mat.location<-which(flea.mat==1, arr.ind=T)
    flea.mat.list<-list(empty.cells=0, flea.mat=flea.mat, 
                        flea.mat.location=flea.mat.location, current.step=0)
    flea.mat.list<-flea.jump(flea.mat.list)
    new.empty.cells.sum<- empty.cells.old + flea.mat.list$empty.cells
    current.iteration<-current.iteration + 1
     if(current.iteration > number.of.repetitions)
      {new.empty.cells<-new.empty.cells.sum/(current.iteration)
       return(new.empty.cells)}
    flea.circus.repetition(flea.mat.list=flea.mat.list, number.of.repetitions, 
                           current.iteration=current.iteration, empty.cells.total=new.empty.cells.sum)}

# Byte-compile the function
flea.circus.repetition<-cmpfun(flea.circus.repetition)
# Utilizing a single CPU
gc()
flea.mat<-matrix(1, nrow=30, ncol=30)
flea.mat.location<-which(flea.mat==1, arr.ind=T)
flea.mat.list<-list(empty.cells=0, 
                    flea.mat=flea.mat,
                    flea.mat.location=flea.mat.location, 
                    current.step=0)

number.of.repetitions<-250

system.time(flea.circus.repetition.out<-flea.circus.repetition(flea.mat.list=flea.mat.list, 
                                                               number.of.repetitions=number.of.repetitions, 
                                                               current.iteration=0,
                                                               empty.cells.total=0))
flea.circus.repetition.out

Function  monte.carlo.rep.func

# For use with parallelization within R 

monte.carlo.rep.func<-function(monte.carlo.reps=1, number.of.repetitions=100)
  {
    number.of.repetitions<-number.of.repetitions
    monte.carlo.reps<-monte.carlo.reps
    replication.mat<-matrix(0, nrow=monte.carlo.reps, ncol=1)
    for(i in 1:monte.carlo.reps)
    {
     gc()
     cat(i,"\n")
     flea.mat<-matrix(1, nrow=30, ncol=30)
     flea.mat.location<-which(flea.mat==1, arr.ind=T)
     flea.mat.list<-list(empty.cells=0, flea.mat=flea.mat, flea.mat.location=flea.mat.location, current.step=0)
     system.time(replication.mat[i, ]<-flea.circus.repetition(flea.mat.list=flea.mat.list,
                                                              number.of.repetitions=number.of.repetitions, 
                                                              current.iteration=0,empty.cells.total=0))
     cat(replication.mat[i, ], "\n")
    }
    return(replication.mat)
  }

# Byte-compile the function
monte.carlo.rep.func<-cmpfun(monte.carlo.rep.func)

Parallelization of  Function monte.carlo.rep.func

# Load libraries

require(foreach)
require(parallel)
require(doMC)

# Register the availability of 8 core CPU 
registerDoMC(8)

gc()
system.time(flea.list.out<-foreach(i=1:8, .verbose=TRUE) %dopar%
       {
          monte.carlo.rep.func(monte.carlo.reps=1, number.of.repetitions=200)
       })

# Extract contents of list structure and take the average of replications across each list
# element
mean(rapply(flea.list.out, as.numeric))

Example Output

   

#  0 replications using 1 CPU 
# 
# 
#    user  system elapsed 
#   1.604   0.012   1.612       ~ 1.6 sec.
# > flea.circus.repetition.out
# [1] 324.5132
# 
# 
# 
# 
# 250 replications using 1 CPU
# 
# 
# 
#    user  system elapsed 
# 358.460   0.080 358.034      ~ 6 min.
# > flea.circus.repetition.out
# [1] 330.9852
# 
#  
#
# 
# 2000 replications using 8 CPUs
#
#
# 
#     user   system  elapsed 
# 4194.288    1.588  381.568   ~ 6.4 min.
# > 
# > mean(rapply(flea.list.out, as.numeric))
# [1] 330.6855