Installation
To install the gridabm library, you need the remotes package
to install it from Github. If you haven’t installed remotes
yet, you can get it from CRAN:
install.packages("remotes")
With the remotes package, you can install
gridabm from Github:
remotes::install_github("socialdynamicshub/gridabm")
First Steps
The gridabm library lets you play around with cellular
automata on 2D grids. You can think of a grid just as a board of
cells.
In the gridabm library, the grid is represented by a
matrix whose entries only take on values from a set of distinct
numeric states. Let’s for example create a matrix with the distinct
states \(0\) and \(1\). Let’s say state \(0\) means that the cell is empty, so we
assign it the color white. State \(1\)
represents “something” on the board, we’ll assign it the color
black.
Then the following matrix…
\[ \begin{pmatrix} 0 & 0 & 0 & 0 & 1\\ 0 & 1 & 0 & 0 & 1\\ 1 & 0 & 0 & 0 & 1\\ 0 & 1 & 1 & 0 & 0\\ 1 & 0 & 0 & 0 & 0 \end{pmatrix} \]
… corresponds to the following plot:
You can construct boards with more distinct states in the same way.
Update Rules
Apart from a grid, a set of distinct states, and an initial state, we need a so-called update rule. In agent-based modeling, this is sometimes referred to as a stepping function. Usually, the update rules are defined at the cell/agent level, meaning that the state of an agent in the next iteration of the automaton is contingent on the states of its neighbors (or a subset of them).
What constitutes the neighborhood of an agent can be defined in different ways, but common neighborhoods are implemented in the package. For instance, the Moore neighborhood is simply the set of all adjacent cells, including diagonally adjacent ones.
The Moore neighborhood of the darker cell is given by the lighter
cells. To get the coordinates of a cell’s Moore neighborhood in the
gridabm package, use the
get_moore_neighborhood function.
get_moore_neighborhood(i = 4, j = 5, dims = c(20, 15), periodic = TRUE)
i and j are the coordinates in the matrix,
axis_size refers to the dimension of the matrix in either
direction1, and the periodic
parameter indicates whether the grid wraps around (i.e., more precisely,
we’re dealing with a torus).
In the gridabm library, update rules are simply
functions that
- take a matrix with numeric states,
- perform the update rule for each cell,
- and return the updated matrix.
The updated matrix should only contain values in the set of distinct states that we defined in the beginning.
Let’s examine this with a very simple example. We use a matrix with distinct states \(0\) and \(1\). Now we define an update rule as follows:
Each cell adopts the state of the cell to its right.
We could implement this rule like this:
Now, let’s create an initial state, run the simulation, and animate it. To change things up, let’s assign the color white to state \(0\) and blue to state \(1\)
left_shift_board <- create_board(2, c(0.7, 0.3), c(20, 20))
left_shift_results <- run_automaton(
left_shift_board,
steps = 100,
stepfunc = left_shift_rule
)
animate_model_run(
left_shift_results,
marker_size = 5,
color_scheme = c(`0` = "white", `1` = "blue")
)
The code above is the general workflow of gridabm.
First, you create a state that you want to start from. It’s a simple
matrix with several distinct states.
You can also create boards more conveniently with the
create_board function. The following code is equivalent to
the one above.
Then, you run the automaton with the stepping function
(stepfunc) that you defined before.
left_shift_results <- run_automaton(
left_shift_board,
steps = 100,
stepfunc = left_shift_rule
)
The results can be animated with the function
animate_model_run.
animate_model_run(
left_shift_results,
marker_size = 5,
color_scheme = c(`0` = "white", `1` = "blue")
)
Color Themes
For the update rules implemented in the gridabm library,
some predefined color schemes for plotting are provided as themeing
functions. For instance, the Schelling model (Schelling 1971) comes with a commonly
used theme of “orange” and “blue” agents.
schelling_board <- create_board(3, c(0.1, 0.45, 0.45), c(20, 20))
plot_state(schelling_board, marker_size = 5, color_scheme = theme_schelling())
Defining your own color schemes is very easy too. As automata states
are simply represented by numeric states, you can just supply a vector
where the color at position i represents state
i. For instance, we could change the colors of the board
above like this:
plot_state(
schelling_board,
marker_size = 5,
color_scheme = c(`0` = "transparent", `1` = "gold2", `2` = "royalblue")
)
The outputs of the plot_state function are
ggplot objects, so you can even theme them further using
the functionality of ggplot2. For
instance, if you would like to plot this board with a dark theme, you
can do the following:
plot_state(
schelling_board,
marker_size = 5,
color_scheme = c(`0` = "transparent", `1` = "gold2", `2` = "royalblue")
) +
theme(panel.background = element_rect(fill = "grey10"))
The same is true for the outputs of animate_model_run.
They are animations created with the gganimate
library, so they can be modified with ggplot2 as
well.
schelling_board %>%
run_automaton(steps = 100, stepfunc = schelling_step, tolerance = 3) %>%
animate_model_run(
marker_size = 5,
color_scheme = c(`0` = "transparent", `1` = "gold2", `2` = "royalblue")
) +
theme(panel.background = element_rect(fill = "grey10"))
Simulating a Forest Fire
Several well-known cellular automata (or agent-based models, depending on how you view it) are implemented in this package. One of them is the so-called forest fire simulation (Drossel and Schwabl 1992).
Let’s first create a board. The color coding is intended to mimic an
actual forest fire: Green cells are trees, orange cells are trees on
fire, dark red cells are burnt trees, and white cells are empty. For the
forest fire, there is a convenience function
create_forest_board because there is a specific way this
model is usually setup. You can see in the plot that the trees on the
left of the board are all on fire and the rest of the trees are
scattered across the board.
create_forest_board(dims = c(20, 20), tree_density = 0.7) %>%
plot_state(marker_size = 5, color_scheme = theme_forest_fire())
In the create_forest_board function, we set a parameter
tree_density. This is just the probability with which a
tree is put on each cell. An example of a sparser forest would look like
this:
create_forest_board(dims = c(20, 20), tree_density = 0.3) %>%
plot_state(marker_size = 5, color_scheme = theme_forest_fire())
… and a more dense one would look like this:
create_forest_board(dims = c(20, 20), tree_density = 0.9) %>%
plot_state(marker_size = 5, color_scheme = theme_forest_fire())
The stepping function for the forest fire is called
forest_fire_step. The rule can be summarized s follows:
Each tree that is on fire sets all its neighboring trees (Von Neumann neighborhood) on fire. Each tree on fire transitions into the state “burnt”.
Note that this is a simplified version of the model described by in Drossel and Schwabl (1992).
You can run and animate the automaton like this:
forest <- create_forest_board(c(20, 20), 0.7)
forest_results <- run_automaton(forest, 100, forest_fire_step)
animate_model_run(forest_results, 5, theme_forest_fire())
As you can see, a wall of fire runs through the forest pretty quickly.
Let’s see, how this looks like in a less dense forest.
forest <- create_forest_board(c(20, 20), 0.5)
forest_results <- run_automaton(forest, 100, forest_fire_step)
animate_model_run(forest_results, 5, theme_forest_fire())
The fire doesn’t run through the entire forest, but stops soon after starting, leaving a significant portion of the forest unburnt.
Now these animations are nice, but we can go beyond merely looking at
the model and analyze the effects of different parameters in-depth. The
run_automaton function returns a data.frame
that contains the results of the model run.
forest <- create_forest_board(c(20, 20), 0.7)
forest_results <- run_automaton(forest, 100, forest_fire_step)
head(forest_results, 20)
| step | cell_id | row | col | state |
|---|---|---|---|---|
| 0 | 1 | 1 | 1 | 0 |
| 0 | 2 | 1 | 2 | 0 |
| 0 | 3 | 1 | 3 | 0 |
| 0 | 4 | 1 | 4 | 0 |
| 0 | 5 | 1 | 5 | 1 |
| 0 | 6 | 1 | 6 | 1 |
| 0 | 7 | 1 | 7 | 0 |
| 0 | 8 | 1 | 8 | 1 |
| 0 | 9 | 1 | 9 | 1 |
| 0 | 10 | 1 | 10 | 1 |
| 0 | 11 | 1 | 11 | 1 |
| 0 | 12 | 1 | 12 | 1 |
| 0 | 13 | 1 | 13 | 0 |
| 0 | 14 | 1 | 14 | 0 |
| 0 | 15 | 1 | 15 | 1 |
| 0 | 16 | 1 | 16 | 1 |
| 0 | 17 | 1 | 17 | 1 |
| 0 | 18 | 1 | 18 | 0 |
| 0 | 19 | 1 | 19 | 1 |
| 0 | 20 | 1 | 20 | 1 |
Let’s look how tree density is related to how quickly the fire runs through the forest and how much of the forest is burnt.
# Run 5 replicates of the model
d <- data.frame()
for (density in seq(0.1, 0.9, 0.1)) {
for (replicate in seq(1, 5, 1)) {
forest <- create_forest_board(c(20, 20), density)
forest_results <- run_automaton(forest, 100, forest_fire_step)
forest_results$density <- density
forest_results$replicate <- replicate
d <- dplyr::bind_rows(d, forest_results)
}
}
# Format the data and plot the results
d %>%
filter(step <= 40) %>%
filter(state %in% c(1, 2, 3)) %>%
mutate(burnt = factor(state == 3, levels = c(FALSE, TRUE))) %>%
group_by(density, replicate, step, burnt, .drop = FALSE) %>%
summarize(count = n()) %>%
ungroup() %>%
group_by(density, replicate, step) %>%
mutate(forest_size = sum(count)) %>%
ungroup() %>%
mutate(relfreq = count / forest_size) %>%
group_by(density, step, burnt) %>%
summarize(mean_relfreq = mean(relfreq)) %>%
ungroup() %>%
ggplot(aes(x = step, y = mean_relfreq, color = burnt)) +
geom_line() +
geom_point(shape = 21, fill = "white") +
scale_y_continuous(breaks = c(0, 0.25, 0.5, 0.75, 1)) +
scale_color_manual(values = c("green", "firebrick")) +
facet_wrap(. ~ density, ncol = 3) +
labs(
color = "Burnt",
x = "Step",
y = "Mean Relative Frequency [5 Replicates]"
) +
theme(
panel.background = element_rect(fill = "white", color = "lightgrey"),
strip.background = element_rect(fill = "white", color = "transparent"),
panel.grid.major.x = element_line(color = "lightgrey"),
panel.grid.major.y = element_line(color = "lightgrey"),
axis.ticks = element_line(color = "lightgrey")
)
This plot displays the average percentage of burnt vs. non-burnt trees in 5 replicates of the model over time. You can clearly observe the general trend that denser forests loose higher percentages of trees. However, in the first three runs (\(density = 0.1\) to \(density = 0.3\)), you can see a reduction in the percentage of trees burnt is reduced. This is due to a subtle error in the data formatting. Can you spot it?
Simulating Urban Segregation
A “classic” agent-based model is the Schelling model of urban segregation (Schelling 1971). Schelling demonstrated that even with seemingly high “tolerance” levels in urban areas2, stark segregation can occur.
We set up the model with agents of two kinds, let’s call them “orange” and “blue” agents.
schelling_board <- create_board(3, c(0.1, 0.45, 0.45), c(20, 20))
plot_state(schelling_board, marker_size = 5, color_scheme = theme_schelling())
The update rule is defined as:
If an agent has more agents of the other kind in its Moore neighborhood than a tolerance parameter \(t\), it relocates to a random free position on the board.
This means that if the tolerance level is set to \(t = 3\), an orange agent will relocate if it has more than \(3\) blue neighbors. Let’s see how this looks like.
schelling_board %>%
run_automaton(steps = 100, stepfunc = schelling_step, tolerance = 3) %>%
animate_model_run(marker_size = 5, color_scheme = theme_schelling())
As you can see, the agents segregate into larger homogeneous areas.
schelling_board %>%
run_automaton(steps = 100, stepfunc = schelling_step, tolerance = 5) %>%
animate_model_run(marker_size = 5, color_scheme = theme_schelling())
If you look at a model run with a tolerance level \(t = 5\), we can still observe substantial segregation, even though the agents have an ostensibly high tolerance level.
There are many things that can be analyzed in this model, but I’ll leave this to you. Here are some primers though.
Manipulate the Distribution of Agents
For instance, what would happen if the orange agents are the \(2:1\) majority? You can control the agent
distribution with the prob argument in the R base
sample function.
schelling_board <- create_board(
n_states = 3,
state_dist = c(0.1, 0.6, 0.3),
dims = c(20, 20)
)
plot_state(schelling_board, marker_size = 5, color_scheme = theme_schelling())
Use More Agent Types
In the previous examples, we used two different kinds of agents, but you can use more agent types as well.
schelling_board <- create_board(4, c(0.1, 0.3, 0.3, 0.3), c(20, 20))
plot_state(
schelling_board,
marker_size = 5,
color_scheme = c(`0` = "transparent", `1` = "gold2", `2` = "royalblue", `3` = "coral")
)
Other Rules Implemented in
gridabm
Conway’s Life
life_board <- create_board(2, c(0.8, 0.2), c(25, 25))
life_results <- run_automaton(life_board, 100, life_step)
animate_model_run(life_results, 3, theme_life())
Brian’s Brain
brians_brain_board <- create_board(3, c(0.5, 0.25, 0.25), c(50, 50))
brians_brain_results <- run_automaton(
brians_brain_board,
200,
brians_brain_step
)
animate_model_run(brians_brain_results, 2, theme_brians_brain_dark()) +
theme(
panel.background = element_rect(fill = "black"),
panel.grid.major = element_line(color = "black"),
panel.grid.minor = element_line(color = "black")
)
Langton’s Ant
langtons_ant_board <- create_ant_board(
antpos = c(36, 36),dims = c(71, 71), antorient = "north"
)
langtons_ant_results <- run_automaton(langtons_ant_board, 200, langtons_ant_step)
animate_model_run(langtons_ant_results, 1, theme_langtons_ant())
After 10,000 steps:
