Animating Schelling’s segregation model

April 29, 2012
By

(This article was first published on R snippets, and kindly contributed to R-bloggers)

Recent blog post on Animations in R inspired me to write a code that generates animations of simulation model. For this task I have chosen Schelling's segregation model.

Having written the code I have found that one year ago a similar code has been proposed. However, the implementation model is different so I thought it is a nice comparison.

Here is the code implementing the model:

# 0 - empty
# 2 - first agent type color
# 4 - second agent type color

# initialize simulation
# size      - square size
# perc.full - percentage of lots to be occupied
init <- function(side, perc.full) {
    size <- floor(side ^ 2 * perc.full / 2)
    state <- matrix(0, side, side)
    occupied <- sample(side ^ 2, 2 * size)
    state[occupied] <- c(2,4)
    return(state)
}

# plot simulation state
# state - simulation state
# i     - simulation iteration
do.plot <- function(state, i) {
    side <- dim(state)[1]
    x <- rep(1:side, side)
    y <- rep(1:side, each = side)
    par(fin=c(4,4), fig=c(0,1,0,1))
    plot(x , y, axes = F, xlab="", ylab="", col = state,
         main = paste("Step", i), pch = 19, cex = 40 / side)
}

# perform one step of simulation
# state     - simulation state
# threshold - percent of required agents of the same color
#             in neighborhood
# radius    - neighborhood radius
sim.step <- function(state, threshold, radius) {
    mod.1 <- function(a, b) { 1 + ((a - 1) %% b) }
    div.1 <- function(a, b) { 1 + ((a - 1) %/% b) }

    unhappy <- rep(NA, length(state))
    side <- dim(state)[1]
    check <- (-radius):(radius)
   
    #find unhappy agents
    for (n in which(state > 0)) {
      x <- div.1(n, side)
      y <- mod.1(n, side)
      x.radius <- mod.1(check + x, side)
      y.radius <- mod.1(check + y, side)
      region <- state[y.radius, x.radius]
      similar <- sum(region == state[n]) - 1
      total <- sum(region > 0) - 1
        unhappy[n] <- (similar < total * threshold)
    }
    vunhappy <- which(unhappy)

    # move unhappy agents
    vunhappy <- vunhappy[sample.int(length(vunhappy))]
    empty <- which(state == 0)
    for (n in vunhappy) {
      move.idx <- sample.int(length(empty), 1)
      state[empty[move.idx]] <- state[n]
      state[n] <- 0
      empty[move.idx] <- n
    }
    return(state)
}

library(animation)

# simple wrapper for animation plotting
go <- function() {
    s <- init(51, 0.75)
    for (i in 1:50) {
        do.plot(s, i)
        last.s <- s
        s <- sim.step(s, 0.6, 1)
        if (identical(last.s, s)) { break }
    }
    for (j in 1:4) {
        do.plot(s, i)
    }
    ani.options(interval = 5 / (i + 2))
}

saveGIF(go())

and a sample animation it generates:


To leave a comment for the author, please follow the link and comment on his blog: R snippets.

R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...



If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Comments are closed.