Schelling's segregation model

In this introductory example we demonstrate Agents.jl's architecture and features through building the following definition of Schelling's segregation model:

• Agents belong to one of two groups (0 or 1).
• The agents live in a two-dimensional grid with a Chebyshev metric. This leads to 8 neighboring positions per position (except at the edges of the grid).
• Each position of the grid can be occupied by at most one agent.
• If an agent has at least 3 neighbors belonging to the same group, then it is happy.
• If an agent is unhappy, it keeps moving to new locations until it is happy.

Schelling's model shows that even small preferences of agents to have neighbors belonging to the same group (e.g. preferring that at least 3/8 of neighbors to be in the same group) could lead to total segregation of neighborhoods.

This model is also available as Models.schelling.

Creating a space

For this example, we will be using a Chebyshev 2D grid, e.g.

using Agents

space = GridSpace((10, 10); periodic = false)

GridSpace with size (10, 10), metric=chebyshev, periodic=false

Agents belonging in this type of space must have a position field that is a NTuple{2, Int}. We ensure this below.

Defining the agent type

mutable struct SchellingAgent <: AbstractAgent
    id::Int             # The identifier number of the agent
    pos::NTuple{2, Int} # The x, y location of the agent on a 2D grid
    mood::Bool          # whether the agent is happy in its position. (true = happy)
    group::Int          # The group of the agent, determines mood as it interacts with neighbors
end

We added two more fields for this model, namely a mood field which will store true for a happy agent and false for an unhappy one, and an group field which stores 0 or 1 representing two groups.

Notice also that we could have taken advantage of the macro @agent (and in fact, this is recommended), and defined the same agent as:

@agent SchellingAgent GridAgent{2} begin
    mood::Bool
    group::Int
end

Creating an ABM

To make our model we follow the instructions of AgentBasedModel. We also want to include a property min_to_be_happy in our model, and so we have:

properties = Dict(:min_to_be_happy => 3)
schelling = ABM(SchellingAgent, space; properties)

AgentBasedModel with 0 agents of type SchellingAgent
 space: GridSpace with size (10, 10), metric=chebyshev, periodic=false
 scheduler: fastest
 properties: Dict(:min_to_be_happy => 3)

Here we used the default scheduler (which is also the fastest one) to create the model. We could instead try to activate the agents according to their property :group, so that all agents of group 1 act first. We would then use the scheduler Schedulers.by_property like so:

schelling2 = ABM(
    SchellingAgent,
    space;
    properties = properties,
    scheduler = Schedulers.by_property(:group),
)

AgentBasedModel with 0 agents of type SchellingAgent
 space: GridSpace with size (10, 10), metric=chebyshev, periodic=false
 scheduler: property
 properties: Dict(:min_to_be_happy => 3)

Notice that Schedulers.by_property accepts an argument and returns a function, which is why we didn't just give Schedulers.by_property to scheduler.

Creating the ABM through a function

Here we put the model instantiation in a function so that it will be easy to recreate the model and change its parameters.

In addition, inside this function, we populate the model with some agents. We also change the scheduler to Schedulers.randomly. Because the function is defined based on keywords, it will be of further use in paramscan below.

using Random # for reproducibility
function initialize(; numagents = 320, griddims = (20, 20), min_to_be_happy = 3, seed = 125)
    space = GridSpace(griddims, periodic = false)
    properties = Dict(:min_to_be_happy => min_to_be_happy)
    rng = Random.MersenneTwister(seed)
    model = ABM(
        SchellingAgent, space;
        properties, rng, scheduler = Schedulers.randomly
    )

    # populate the model with agents, adding equal amount of the two types of agents
    # at random positions in the model
    for n in 1:numagents
        agent = SchellingAgent(n, (1, 1), false, n < numagents / 2 ? 1 : 2)
        add_agent_single!(agent, model)
    end
    return model
end

initialize (generic function with 1 method)

Notice that the position that an agent is initialized does not matter in this example. This is because we use add_agent_single!, which places the agent in a random, empty location on the grid, thus updating its position.

Defining a step function

Finally, we define a step function to determine what happens to an agent when activated.

function agent_step!(agent, model)
    minhappy = model.min_to_be_happy
    count_neighbors_same_group = 0
    # For each neighbor, get group and compare to current agent's group
    # and increment count_neighbors_same_group as appropriately.
    # Here `nearby_agents` (with default arguments) will provide an iterator
    # over the nearby agents one grid point away, which are at most 8.
    for neighbor in nearby_agents(agent, model)
        if agent.group == neighbor.group
            count_neighbors_same_group += 1
        end
    end
    # After counting the neighbors, decide whether or not to move the agent.
    # If count_neighbors_same_group is at least the min_to_be_happy, set the
    # mood to true. Otherwise, move the agent to a random position.
    if count_neighbors_same_group ≥ minhappy
        agent.mood = true
    else
        move_agent_single!(agent, model)
    end
    return
end

agent_step! (generic function with 1 method)

For the purpose of this implementation of Schelling's segregation model, we only need an agent step function.

When defining agent_step!, we used some of the built-in functions of Agents.jl, such as nearby_positions that returns the neighboring position on which the agent resides, ids_in_position that returns the IDs of the agents on a given position, and move_agent_single! which moves agents to random empty position on the grid. A full list of built-in functions and their explanations are available in the API page.

Stepping the model

Let's initialize the model with 370 agents on a 20 by 20 grid.

model = initialize()

AgentBasedModel with 320 agents of type SchellingAgent
 space: GridSpace with size (20, 20), metric=chebyshev, periodic=false
 scheduler: randomly
 properties: Dict(:min_to_be_happy => 3)

We can advance the model one step

step!(model, agent_step!)

Or for three steps

step!(model, agent_step!, 3)

Running the model and collecting data

We can use the run! function with keywords to run the model for multiple steps and collect values of our desired fields from every agent and put these data in a DataFrame object. We define vector of Symbols for the agent fields that we want to collect as data

adata = [:pos, :mood, :group]

model = initialize()
data, _ = run!(model, agent_step!, 5; adata)
data[1:10, :] # print only a few rows

10 rows × 5 columns

We could also use functions in adata, for example we can define

x(agent) = agent.pos[1]
model = initialize()
adata = [x, :mood, :group]
data, _ = run!(model, agent_step!, 5; adata)
data[1:10, :]

10 rows × 5 columns

With the above adata vector, we collected all agent's data. We can instead collect aggregated data for the agents. For example, let's only get the number of happy individuals, and the average of the "x" (not very interesting, but anyway!)

using Statistics: mean
model = initialize();
adata = [(:mood, sum), (x, mean)]
data, _ = run!(model, agent_step!, 5; adata)
data

6 rows × 3 columns

Other examples in the documentation are more realistic, with more meaningful collected data. Don't forget to use the function dataname to access the columns of the resulting dataframe by name.

Visualizing the data

We can use the abm_plot function to plot the distribution of agents on a 2D grid at every generation, via the InteractiveDynamics.jl package and the Makie.jl plotting ecosystem.

Let's color the two groups orange and blue and make one a square and the other a circle.

using InteractiveDynamics
using CairoMakie # choosing a plotting backend

groupcolor(a) = a.group == 1 ? :blue : :orange
groupmarker(a) = a.group == 1 ? :circle : :rect
figure, _ = abm_plot(model; ac = groupcolor, am = groupmarker, as = 10)
figure # returning the figure displays it

Animating the evolution

The function abm_video can be used to save an animation of the ABM into a video. You could of course also explicitly use abm_plot in a record loop for finer control over additional plot elements.

model = initialize();
abm_video(
    "schelling.mp4", model, agent_step!;
    ac = groupcolor, am = groupmarker, as = 10,
    framerate = 4, frames = 20,
    title = "Schelling's segregation model"
)

Launching the interactive application

Given the definitions we have already created for a normally plotting or animating the ABM it is almost trivial to launch an interactive application for it, through the function abm_data_exploration.

We define a dictionary that maps some model-level parameters to a range of potential values, so that we can interactively change them.

parange = Dict(:min_to_be_happy => 0:8)

Dict{Symbol, UnitRange{Int64}} with 1 entry:
  :min_to_be_happy => 0:8

We also define the data we want to collect and interactively explore, and also some labels for them, for shorter names (since the defaults can get large)

adata = [(:mood, sum), (x, mean)]
alabels = ["happy", "avg. x"]

model = initialize(; numagents = 300) # fresh model, noone happy

AgentBasedModel with 300 agents of type SchellingAgent
 space: GridSpace with size (20, 20), metric=chebyshev, periodic=false
 scheduler: randomly
 properties: Dict(:min_to_be_happy => 3)

figure, adf, mdf = abm_data_exploration(
    model, agent_step!, dummystep, parange;
    ac = groupcolor, am = groupmarker, as = 10,
    adata, alabels
)

Ensembles and distributed computing

We can run ensemble simulations and collect the output of every member in a single DataFrame. To that end we use the ensemblerun! function. The function accepts a Vector of ABMs, each (typically) initialized with a different seed and/or agent distribution. For example we can do

models = [initialize(seed = x) for x in rand(UInt8, 3)];

and then

adf, = ensemblerun!(models, agent_step!, dummystep, 5; adata)
adf[(end - 10):end, :]

11 rows × 4 columns

It is possible to run the ensemble in parallel. For that, we should start julia with julia -p n where n is the number of processing cores. Alternatively, we can define the number of cores from within a Julia session:

using Distributed
addprocs(4)

For distributed computing to work, all definitions must be preceded with @everywhere, e.g.

using Distributed
@everywhere using Agents
@everywhere mutable struct SchellingAgent ...
@everywhere agent_step!(...) = ...

Then we can tell the ensemblerun! function to run the ensemble in parallel using the keyword parallel = true:

adf, = ensemblerun!(models, agent_step!, dummystep, 5; adata, parallel = true)

Scanning parameter ranges

We often are interested in the effect of different parameters on the behavior of an agent-based model. Agents.jl provides the function paramscan to automatically explore the effect of different parameter values.

We have already defined our model initialization function as initialize. We now also define a processing function, that returns the percentage of happy agents:

happyperc(moods) = count(moods) / length(moods)
adata = [(:mood, happyperc)]

parameters = Dict(
    :min_to_be_happy => collect(2:5), # expanded
    :numagents => [200, 300],         # expanded
    :griddims => (20, 20),            # not Vector = not expanded
)

adf, _ = paramscan(parameters, initialize; adata, agent_step!, n = 3)
adf

32 rows × 4 columns

We nicely see that the larger :min_to_be_happy is, the slower the convergence to "total happiness".