Site icon R-bloggers

Animating Schelling’s segregation model

[This article was first published on R snippets, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
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< o:p>
# 2 – first agent type color< o:p>
# 4 – second agent type color< o:p>

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

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

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

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

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

library(animation)< o:p>

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

saveGIF(go())< o:p>

and a sample animation it generates:


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

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.