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**())**
Hi! Thank you for your code. Can you tell me how you define radius? If you have a radius of 1, is that a radius of 1 agent? And does that mean each agent has 4 neighbors (left, right, up, down) in their radius of 1 or 8 neighbors (left, right, up, down, plus 4 corner agents)?

ReplyDelete8 neighbors (so called Moore neighborhood)

DeleteMaybe I am a bit stupid but I cannot run it in R, I get many errors, the first one being:

ReplyDeleteError: unexpected input in:

" par(fin=c(4,4), fig=c(0,1,0,1))

plot(x , y, axes = F, xlab=“"

I have just tested copying the code from the website and pasting it to R and it just works.

DeleteFrom what you have pasted it seems that you perform copy-paste incorrectly. Observe that the first quotation mark is different than the second after "xlab="