Model Thinking- Schelling Process Functions

SchellingProcessFunctions.R

# these functions were already presented in gist 1885003, but are copied here
# so that they can be sourced directly from the R console.
# to do so, click on the 'raw' button, get rid of the 's' in https, and change the .R to .r, like this:
# source("http://raw.github.com/gist/1910183/cbc854b9f1de4495db6fbd7eba43c337a1a8922d/SchellingProcessFunctions.r")

queenNeighbors <- function(ind, M, n){
	inds <- ind+outer(c(-1, 0, 1), c(-n, 0, n), "+")
	M[inds[inds > 0 & inds <= length(M) & inds != ind]]
}

# determine if an individual of position index 'ind' will move
MoveYes <- function(ind, M, n, threshold = threshold){
	status <- M[ind]
	neighbors <- queenNeighbors(ind = ind, M = M, n = n)
	psame <- sum(neighbors == status, na.rm = TRUE) / length(neighbors[!is.na(neighbors)])
	if (!(is.null(psame)|is.na(psame))){
		if (psame < threshold) return(ind)
	}
}

# determine all who will move
Movers <- function(M, n, threshold = threshold){
	unlist(sapply((1:length(M))[!is.na(M)], MoveYes, M = M, threshold = threshold, n = 10))
}

# carries out a single iteration:
# 1) determine who will move
# 2) what are their statuses
# 3) where do they go (random assignment)
# 4) empty previous occupied spaces
# 5) if at some point no one wants to move, pass on a message to next function to
#    it knows to break the cycle
MoveCycle <- function(M, n = 10, threshold = threshold, nvacants = nvacants){
	movers.all.ind <- Movers(M, n = n, threshold = threshold)
	if (length(movers.all.ind) > 0){
		if (length(movers.all.ind) > nvacants){
			movers.all.ind <- sample(movers.all.ind, nvacants, replace = FALSE)
		}
		movers.all.statuses <- M[movers.all.ind]
		# random reassignment of places movers.all.statuses to available spots:
		M[sample((1:length(M))[is.na(M)], length(movers.all.ind),replace=FALSE)] <- movers.all.statuses
		M[movers.all.ind] <- NA
		return(M)
	} else {
		return(list("breaktheloop", M = M))
	}
}
#    end internal functions ------------------------------------------------------------#

# main function ------------------------------------------------------------------------#
# given starting matrix M0, a 'sameness' threshold and a maximum numer of iterations, 
# walk through the moving process. If at some point everyone is satisfied, then the loop
# breaks
SchellingProcess <- function(M0, threshold = .3, maxit = 100){
	nvacants <- sum(is.na(M0))
	M <- list()
	M[[1]] <- M0
	n <- ncol(M)
	for (i in 1:maxit){
		M[[i+1]] <- MoveCycle(M[[i]], n = n, threshold = threshold, nvacants = nvacants)
		if (!is.matrix(M[[i + 1]])){
			break
		}
	}
	out <- list(M = M, stable = !(i == maxit), threshold = threshold)
	class(out) <- "Schelling"
	return(out)
}
# end main function --------------------------------------------------------------------#

# define plot method -------------------------------------------------------------------#
# plot() method for output from SchellingProcess(), which is of class "Schelling"
plot.Schelling <- function(M, frame.pause = .1,col=c("darkblue", "green3")){
	iterations <- length(M$M)-1
	for (i in 1:iterations){
		Sys.sleep(frame.pause)
		image(t(M$M[[i]]), useRaster = TRUE, col = col, axes = FALSE)
	}
}
# end plot method ----------------------------------------------------------------------#

# function to make starting population -------------------------------------------------#
# just an easy way to get a starting matrix randomly filled in
PopulateMatrix <- function(nrow = 10, ncol = 10, prob = c(.3,.7), propvacant = .2){
	l <- nrow*ncol
	samplesize <- round(l * (1 - propvacant))
	M <- matrix(NA, nrow = nrow, ncol = ncol)
	M[sample(1:prod(dim(M)), samplesize, replace = FALSE)] <- sample(c(TRUE, FALSE), samplesize, prob = c(.3, .7), replace = TRUE)
	return(M)
}
# --------------------------------------------------------------------------------------#