I want to play a game ... Conway's Game of Life! For explanation, see http://en.wikipedia.org/wiki/Conway's_Game_of_Life

GameOfLife.R

#' Execute a round of the Game of Life
#' 
#' I want to play a game. Specifically, Conway's Game of Life. 
#' 
#' @param x a matrix populated with 0s and 1s.
#' @param birth a vector indicating the birthing rule. Defaults to 3.
#' @param stay a vector indicating the stay alive rule Defaults to c(2,3).
#' @return a matrix representing an updated input matrix according to the rules
#' @references \url{http://en.wikipedia.org/wiki/Conway's_Game_of_Life}
#' @author Bryan Goodrich
tick <- function(x, birth = 3, stay = c(2, 3)) {
    stopifnot(all(x %in% c(0,1)) || !is.null(birth) || !is.null(stay))

    LIVE    <- 1
    DEAD    <- 0
    ROWSIZE <- dim(x)[1]
    COLSIZE <- dim(x)[2]
    
    # For Corner Wrapping
    fixer <- function(n) {
        function(x) {
            b <- (x-1):(x+1)
            b <- ifelse (b == 0, n, b)
            b <- ifelse (b > n, 1, b)
            return(b)
        }
    }
    
    rule <- function(r) {
        function(x) ifelse (any(x %in% r), LIVE, DEAD)
    }
    
    Rband <- fixer(ROWSIZE)
    Cband <- fixer(COLSIZE)
    B     <- rule(birth)
    SA    <- rule(stay)
    
    
    newx <- x 
    for (i in seq(ROWSIZE)) {
        for (j in seq(COLSIZE)) {
            window    <- sum(x[Rband(i), Cband(j)])
            newx[i,j] <- ifelse (x[i,j] == LIVE, SA(window-1), B(window))
        }
    }
    
    newx
}



# Print function. 
# -- Probably more appropriate to have tick return a classed matrix with a print/plot version here.
display_tick <- function(m, xaxt = 'n', yaxt = 'n', 
                         col = c("blanchedalmond", "darkgreen"), ...) {
    image(t(m), xaxt = xaxt, yaxt = yaxt, col = col, ...)
}


# Gosper's Glider Gun Example
glider <- matrix(inverse.rle(structure(list(lengths = c(1526L, 2L, 69L, 2L, 637L, 3L, 67L, 
1L, 3L, 1L, 65L, 1L, 5L, 1L, 64L, 1L, 5L, 1L, 67L, 1L, 68L, 1L, 
3L, 1L, 67L, 3L, 69L, 1L, 209L, 3L, 68L, 3L, 67L, 1L, 3L, 1L, 
136L, 2L, 3L, 2L, 705L, 2L, 69L, 2L, 1527L), values = c(0, 1, 
0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 
1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 
0)), .Names = c("lengths", "values"), class = "rle")), ncol = 78)

m <- glider
for (n in 1:500) {
    display_tick(m)
    m <- tick(m)
    Sys.sleep(0.01)  # May not be needed for you to see it.
}



# B36/S23 "High Life" Alternative
# Set your own seed and look for Replicators!
# http://www.youtube.com/watch?v=Mw-YPFlPv2U
set.seed(666)
N <- 100
p <- 0.3

m <- sample(c(0, 1), N*N, replace = TRUE, prob = c(1-p, p))
dim(m) <- c(N, N)
for (n in 1:500) {
    display_tick(m)
    m <- tick(m, birth = c(3, 6))  
}

Thanks for the insight! I'm not very experienced with filter, so I'll have to toy around with exactly what is being done in those 3 lines to fully appreciate it. I suspect a more vectorized solution to what I've done will be faster, which is a great limitation in my code. My initial plan was to do something in Rcpp to handle the "heavy lifting." I still plan to do that, but the code was designed with that sort of logic in mind. I never intended to do a pure R implementation other than to see how this cellular automata worked.

Hadley Wickham had a good idea that it would be nice to have an import mechanism for reading in a file (string) of text that indicates an initial matrix. That shouldn't be too hard to implement, but I also think an XLSX import will be more useful because then you can manually set up a spreadsheet with the initial state more easily and import it into the game.

One thing about your final statement is that it depends on that exact specification of the birth and stay rule. They can come in many forms, whether sequential (e.g.: 1, 2, 3) or not (e.g.: 1, 5, 8) and the birth rule need not be singular--here an all(x %in% birth) type expression will work.

I also realize that since I'm dealing with integers, I should be using the 0L and 1L instead of 0 and 1.

Again, thank you for the contribution!

Hi Bryan,

filter() takes a (usually long) vector and the filter (usually short) and performs the inner product, with the filter sliding along the vector. For example, if the filter is (1,1,1), then the filtering operation is just adding each three consecutive vector entries; if the filter is (1,0,1), then the sum skips the middle element. The R function was created for time series, and can work with multivariate time series objects, which in some sense are just like matrices (the filtering operation is performed on each column). The operation can be made "circular," which nicely takes care of the edge wrapping.

Linear algebra operations in R (like filtering) are usually coded in C, and make use of optimized linear algebra libraries, so they tend to be much faster than explicit R loops. Now, since you are going to translate your code into C++, those loops will probably be quite efficient too; still, taking advantage of optimized libraries can give an extra edge.

About the final statement, indeed it does not work in full generality, but the modification you suggest should work. ifelse() is quite nice because it is vectorized, and even more, preserves matrix structure.

This is fun! I hadn't played with life in a long time. If one sets stay=6, it tends to produce fingerprint-like patterns!

Cheers,

Omar.

I've used filter on univariate time series for smoothing, but this does make me appreciate it in a new fashion. So on a matrix, it just treats each column as a separate series to be filtered. I used image on each of the intermediary c1, c2, ..., to see how it changed the image of the provided matrix. Very cool! I'm still planning to work on a Coursera class that did a bit of signal processing (MATLAB). That should help. Still, without working out a small sample myself (on my todo list), I've no intuition with if that slight of hand is producing the same result for each of the neighborhoods. I guess I could run a check, of course! It is way faster, though. I like this!

As for filter returning a ts matrix, you can just unclass it and it'll strip the attributes and classes associated with it, returning a vector with the appropriate dimensions (i.e., a raw matrix).

I still plan to do my C++ implementation, because I want to do something more graph and object oriented that'll be more complex (and make room for more complex games!), but this is definitely a great improvement. Thank you so much!