Agent Based Models - Scientific Computing with Python

An agent-based model is a way of modeling some sort of phenomenon using discrete “agents” which interact with other agents, sometimes in very complex ways. In case this is too abstract, you can think of agents as simulated people, and you’re trying to model how they interact with each other.

Example: Spreading an Idea¶

Let’s say someone has a new idea, and that it is a really good idea. So good in fact that once someone hears this new idea they can’t stop talking about it. We can model a person, or agent as a Python class with a single attribute:

class Person():
    """
    An agent representing a person.
    
    By default, a person is not enlightened.  They can become enlightened using the englighten method.
    """
    
    def __init__(self):
        self.enlightened = False # if enlightened = True, the person has heard idea

        
    def is_enlightened(self):
        """
        returns true if the person has heard the idea
        """
        return self.enlightened

    
    def enlighten(self):
        """
        once the person has heard the idea, they become enlightened
        """
        self.enlightened = True
        
        
a = Person()
print(a.is_enlightened()) # the person is not enlightened by default
a.enlighten() # the person becomes enlightened
print(a.is_enlightened()) # the person is now enlightened

False
True

import numpy as np
from numpy.random import randint, rand
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation

from tqdm import tqdm

Now, what we want to do is simulate how this idea spreads through a population. Say that there are N people in the population, and each day each person talks to k random other people. If someone talks to someone else who is enlightened, then the new person becomes enlightened with probability p. We’ll start with one enlightened person and keep track of how the idea spreads

N = 1_000 # population size
k = 5     # number of people each person talks to
p = 0.1   # probability that someone is enlightened

pop = [Person() for i in range(N)] # our population
pop[0].enlighten()

now let’s write a function that will simulate all interactions in a day

from numpy.random import randint, rand

for i in range(N):
    if pop[i].is_enlightened():
        # person i enlightens all their contacts
        contacts = randint(N, size=k)
        for j in contacts:
            if rand() < p:
                pop[j].enlighten()

we want to count how many people are enlightened at the end of the day

def count_enlightened(pop):
    return sum(p.is_enlightened() for p in pop)

count_enlightened(pop)

2

now let’s see how the enlightened population changes over time

T = 100 # number of days to simulate
counts = [count_enlightened(pop)]
for t in range(T):
    # update the population
    for i in range(N):
        if pop[i].is_enlightened():
            # person i enlightens all their contacts
            contacts = randint(N, size=k)
            for j in contacts:
                if rand() < p:
                    pop[j].enlighten()
                    
    # add to our counts
    counts.append(count_enlightened(pop))

import matplotlib.pyplot as plt

plt.plot(counts)
plt.show()

Exercises¶

How does changing N, p and k affect the simulation?
How would you change the simulation model to say that there is a probablity q that an enlightened person will become un-enlightened every day?

## Your code here

T = 100 # number of days to simulate
counts = [count_enlightened(pop)]
for t in range(T):
    # update the population
    for i in range(N):
        if pop[i].is_enlightened():
            # person i enlightens all their contacts
            contacts = randint(N, size=k)
            for j in contacts:
                if rand() < p:
                    pop[j].enlighten()
            if rand() < q:
                #pop[j].unenelighten() # sets enligtened to false
                pop[j].enlightened = False
                    
                    
    # add to our counts
    counts.append(count_enlightened(pop))

Phase Diagrams¶

In our example, there are a couple of parameters we can tweak. We’ll focus on p and k. One question we might ask is how the outcome changes as we change these parameters. We’ll count the total number of people who are enlightened on day 20.

def run_simulation(p, k, N=1_000, T=20):
    """
    return the number of people enlightened at time T
    """
    
    pop = [Person() for i in range(N)] # our population
    pop[0].enlighten()
    
    for t in range(T):
        # update the population
        for i in range(N):
            if pop[i].is_enlightened():
                # person i enlightens all their contacts
                contacts = randint(N, size=k)
                for j in contacts:
                    if rand() < p:
                        pop[j].enlighten()
                    
    return count_enlightened(pop)

run_simulation(0.1, 5)

974

Now, let’s run simulations on a range of values for k and p

import numpy as np
from tqdm import tqdm

ks = np.arange(1, 11, dtype=np.int64)
ps = np.logspace(-2,0, 10)

cts = np.zeros((len(ks), len(ps)))
for i, k in tqdm(enumerate(ks)):
    for j, p in enumerate(ps):
        cts[i,j] = run_simulation(p, k)

10it [00:15,  1.60s/it]

plt.figure(figsize=(10,5))
plt.imshow(cts, extent=[np.min(ps), np.max(ps), np.max(ks), np.min(ks)])
plt.colorbar()
# plt.axis('square')
plt.xscale('log')
plt.xlabel('p')
plt.ylabel('k')

plt.show()

from the above, we see that if k and p are large enough that the whole population is enlightened on day 20, and if they are too small then very few are enlightened. There is a phase transition between the blue and yellow regions above, meaning the choice of parameters can cause abrupt changes in the measured behavior of the system

The Game of Life¶

You don’t necessarily need to use Python objects to encode everything about an agent-based model. For instance, you might just use an array to store the state of objects in simple simulations, even when behavior is nonlinear.

We’ll take a look at Conway’s Game of Life implemented using a numpy array. The game of life takes place on a m by n grid, and each element of the grid either contains life or it doesn’t. The game models how life might spread over the area of a continent. The cells are updated at each time step using the following rules:

Any live cell with 0 or 1 live neighbors dies (under-population)
Any live cell with 2 or 3 neighbors survives until the next time step (sustainable population)
Any live cell with more than 3 live neighbors dies (overpopulation)
Any dead cell with exactly 3 neighbors becomes alive (reproduction)

Neighbors are either horizontally, vertically, or diagonally adjacent. We’ll encode whether a cell is alive or dead using 1 and 0 respectively (or True/False)

m = 100
n = 100

S = np.zeros((m,n), dtype=bool) # grid

def neighbors(i, j, m, n):
    inbrs = [-1, 0, 1]
    if i == 0:
        inbrs = [0, 1]
    if i == m-1:
        inbrs = [-1, 0]
    jnbrs = [-1, 0, 1]
    if j == 0:
        jnbrs = [0, 1]
    if j == n-1:
        jnbrs = [-1, 0]
        
    for delta_i in inbrs:
        for delta_j in jnbrs:
            if delta_i == delta_j == 0:
                continue
            yield i + delta_i, j + delta_j

def count_alive_neighbors(S):
    m, n = S.shape
    cts = np.zeros(S.shape, dtype=np.int64)
    for i in range(m):
        for j in range(n):
            for i2, j2 in neighbors(i, j, m, n):
                cts[i,j] = cts[i,j] + S[i2, j2]
                
    return cts

count_alive_neighbors(S)

array([[0, 0, 0, ..., 0, 0, 0],
       [0, 0, 0, ..., 0, 0, 0],
       [0, 0, 0, ..., 0, 0, 0],
       ...,
       [0, 0, 0, ..., 0, 0, 0],
       [0, 0, 0, ..., 0, 0, 0],
       [0, 0, 0, ..., 0, 0, 0]])

What makes the game of life interesting is that you can get patterns that propagate over your grid

m, n = 50, 50
np.random.seed(0)
S = np.random.rand(m, n) < 0.3
plt.figure(figsize=(10,10))
plt.imshow(S)
plt.show()

The below cell produces an animated GIF of the game of life. See here for an example from matplotlib.

from matplotlib.animation import FuncAnimation

fig = plt.figure(figsize=(5,5))
fig.set_tight_layout(True)

# Plot an image that persists
im = plt.imshow(S, animated=True)
plt.axis('off') # turn off ticks

def update(*args):

    global S
    
    # Update image to display next step
    cts = count_alive_neighbors(S)
    # Game of life update
    S = np.logical_or(
        np.logical_and(cts == 2, S),
        cts == 3
    )
    im.set_array(S)
    return im,

anim = FuncAnimation(fig, update, frames=100, interval=200, blit=True)
anim.save('life.gif', dpi=80, writer='imagemagick')
;

''

Here’s what the GIF looks like: