Skip to main content

Beyond disease

🧠10% 🤖90%

Starsim is built for infectious diseases, but its machinery for agents is general: a population of agents, each with properties kept as arrays over all of them; modules that change those properties at each step; results recorded over time; and random numbers drawn so that runs can be repeated and compared. A module written as code can hold any property of people and any rule for changing it, so Starsim can run agent-based models with no disease at all. This page shows two classic ones: a model of wealth, from economics, and a model of flocking birds, from computer graphics.

Not every model of something other than an infection needs code. States and transitions describe any process in which people move between categories at given rates (smoking, use of contraception, stages of a chronic disease, employment), and a declared model of such a process runs under every solver, with the same interventions as a disease: see Diseases. Some models track things that are not people at all, such as bacteria in water or a population of mosquitoes, as states with people=False in the same model as the people (see Cholera with a water reservoir).

Wealth

This is the “wealth model” of the Mesa tutorial: each person starts with one unit of wealth, and gives one unit at a time to someone chosen at random. Though the rule is fair, wealth ends up very unequal. Here each person who has any gives with a chance of 1 in 10 at each step, so that inequality builds up over steps rather than all in the first.

Parameter Value Meaning
People 1,000
Starting wealth 1 unit The same for everyone
Giving 1 unit With a chance of 0.1 a step, from each person who has any, to someone chosen at random
Steps 100

Each person’s wealth is a state, and who gives and to whom are drawn from Starsim distributions, so that runs can be repeated and compared (see Random numbers). The Gini coefficient (0 if everyone has the same, near 1 if one person has everything) is a result:

import numpy as np
import matplotlib.pyplot as plt
import starsim as ss

class Wealth(ss.Module):
    """ Each person starts with 1 unit of wealth; each step, a tenth of those with any give 1 unit to someone chosen at random """

    def __init__(self, **kwargs):
        super().__init__()
        self.define_pars(
            p_give = ss.bernoulli(p=0.1), # The chance that a person who has any wealth gives a unit in a step
        )
        self.update_pars(**kwargs)
        self.define_states(
            ss.FloatArr('wealth', default=1.0, label='Wealth'),
        )
        self.recipient = ss.randint(low=0, high=lambda self, sim, uids: len(sim.people)) # Whom each giver gives to: a position, the same as a UID only since nobody is born or dies here
        return

    def init_results(self):
        super().init_results()
        self.define_results(
            ss.Result('gini', label='Gini coefficient', kind='ratio'),
        )
        return

    def step(self):
        givers = self.pars.p_give.filter((self.wealth > 0).uids)
        recipients = self.recipient.rvs(givers)
        self.wealth[givers] = self.wealth[givers] - 1
        self.wealth[:] = self.wealth + np.bincount(recipients, minlength=len(self.wealth)) # Each recipient gets 1 unit for each giver who chose them
        return

    def record_states(self): # The state at the start of each step
        super().record_states()
        w = np.sort(self.wealth.values)
        n = len(w)
        self.results.gini[self.ti] = np.sum((2*np.arange(1, n+1) - n - 1)*w)/(n*w.sum())
        return
sim = ss.Sim(Wealth(), n_agents=1000, stop=100) # 100 steps
sim.run()
sim.results.wealth.gini.plot()
Sim: wealth, as 1,000 agents, steps 0 to 100
  Done in 0.024 s

Inequality rises over the first 10 or so steps, and then stays high, with a Gini coefficient of about 0.65, close to that of the most unequal countries. The final distribution has many people with nothing and a few with several units:

wealth = sim.people.wealth.wealth
plt.figure()
plt.hist(wealth.values, bins=np.arange(wealth.max() + 2) - 0.5)
plt.xlabel('Wealth')
plt.ylabel('Number of people')
plt.show()

Compared with Mesa, where each agent is an object with its own step(), Starsim keeps each property as one array over all people, so a rule is applied to everyone at once, as operations on arrays.

Flocking

Boids are the classic model of flocking, in which each bird follows three simple rules about the birds it can see: alignment (steer towards their average heading), cohesion (steer towards their center), and separation (steer away from any that are too close). No bird leads, and nobody tells the flock where to go, but flocks form and fly together.

Parameter Value Meaning
Boids 200
Sight 0.1 How far a boid can see, in a world that is a unit square
Too close 0.03 The distance within which a boid steers away
Alignment, cohesion, separation 0.1, 0.005, 3 × 10⁻⁵ How strongly each rule steers
Speed 0.005 to 0.01 The minimum and maximum distance moved in a step
Steps 200

In Starsim, each boid’s position (x, y) and velocity (vx, vy) are states, starting at random. The world is a unit square whose edges wrap around, so a boid leaving on the right comes back on the left. At each step, the module works out the offset from every boid to every other (the short way round), applies the three rules to all the boids at once, keeps each boid’s speed between a minimum and a maximum, and moves them. How well the flock is aligned is a result: for each boid, the length of the average heading of those it can see (itself included), averaged over all boids, which is 1 when each moves with its neighbors and close to 0 when their headings are random:

class Boids(ss.Module):
    """ Flocking: each boid steers towards the heading and center of the boids it can see, and away from those too close """

    def __init__(self, **kwargs):
        super().__init__()
        self.define_pars(
            radius     = 0.1,    # How far a boid can see (the world is a unit square)
            too_close  = 0.03,   # How close is too close
            alignment  = 0.1,    # How strongly a boid steers towards its neighbors' heading
            cohesion   = 0.005,  # ... towards its neighbors' center
            separation = 3e-5,   # ... away from those too close
            speed      = [0.005, 0.01], # The minimum and maximum distance moved in a step
        )
        self.update_pars(**kwargs)
        self.define_states(
            ss.FloatArr('x',  default=ss.random(), label='Position (x)'),
            ss.FloatArr('y',  default=ss.random(), label='Position (y)'),
            ss.FloatArr('vx', default=ss.uniform(-0.01, 0.01), label='Velocity (x)'),
            ss.FloatArr('vy', default=ss.uniform(-0.01, 0.01), label='Velocity (y)'),
        )
        return

    def init_results(self):
        super().init_results()
        self.define_results(ss.Result('aligned', label='Alignment with neighbors', kind='ratio'))
        return

    def offsets(self):
        """ The offsets from each boid to every other, the short way round the edges, and the distances """
        x, y = self.x.values, self.y.values
        dx = x[None, :] - x[:, None]
        dy = y[None, :] - y[:, None]
        dx -= np.round(dx)
        dy -= np.round(dy)
        return dx, dy, np.hypot(dx, dy)

    def step(self):
        p = self.pars
        vx, vy = self.vx.values, self.vy.values
        dx, dy, dist = self.offsets()
        np.fill_diagonal(dist, np.inf) # A boid is not its own neighbor
        near = dist < p.radius
        n_near = np.maximum(near.sum(axis=1), 1)
        close = dist < p.too_close

        # The three rules, for all the boids at once
        steer_x = p.alignment*(near @ vx / n_near - vx) + p.cohesion*(near*dx).sum(axis=1)/n_near - p.separation*(close*dx/dist**2).sum(axis=1)
        steer_y = p.alignment*(near @ vy / n_near - vy) + p.cohesion*(near*dy).sum(axis=1)/n_near - p.separation*(close*dy/dist**2).sum(axis=1)
        alone = ~near.any(axis=1) # A boid that sees nobody flies straight on
        vx = vx + np.where(alone, 0, steer_x)
        vy = vy + np.where(alone, 0, steer_y)

        # Keep the speed within limits, then move, wrapping around the edges
        speed = np.hypot(vx, vy)
        factor = np.clip(speed, *p.speed)/speed
        self.vx[:] = vx*factor
        self.vy[:] = vy*factor
        self.x[:] = (self.x + self.vx) % 1
        self.y[:] = (self.y + self.vy) % 1
        return

    def record_states(self): # The state at the start of each step
        super().record_states()
        speed = np.hypot(self.vx.values, self.vy.values)
        hx, hy = self.vx.values/speed, self.vy.values/speed # Each boid's heading
        _, _, dist = self.offsets()
        sees = dist < self.pars.radius # Including itself
        self.results.aligned[self.ti] = (np.hypot(sees @ hx, sees @ hy)/sees.sum(axis=1)).mean()
        return

The rule of alignment is what makes a flock: without it, boids still gather and keep their distance, but each flies its own way. Here 200 boids fly for 200 steps, with and without it:

scenarios = {
    'Flocking': ss.Sim(Boids(), n_agents=200, stop=200),
    'No alignment': ss.Sim(Boids(alignment=0), n_agents=200, stop=200),
}
msim = ss.parallel(scenarios, workers=1)
msim.plot('boids.aligned')
MultiSim: 2 runs of 2 scenarios (Flocking, No alignment)
  Done in 0.55 s

With alignment, most boids are moving with their neighbors within about 50 steps, and nearly all within 100; without it, their headings stay about as random as at the start.

The positions and velocities are states, so the flock can be drawn at any point in a run. Here the sim runs in stages (run(until=...) runs up to a given step, and the next call carries on from there), keeping the boids’ positions and headings every 3 steps, and then plays them back as an animation, with each boid as an arrow along its heading (the plotting code is folded; click to show it):

The code for the animation 
from matplotlib import animation
from IPython.display import HTML

sim = ss.Sim(Boids(), n_agents=200, stop=200)
sim.init()
boids = sim.custom.boids # The module is given as `custom`: any module that is not a disease, network, demographic module, intervention, analyzer, or connector
frames = []
for step in range(0, 201, 3):
    sim.run(until=step)
    frames.append([state.values.copy() for state in [boids.x, boids.y, boids.vx, boids.vy]])

fig, ax = plt.subplots(figsize=(4, 4), dpi=72)
ax.set(xlim=(0, 1), ylim=(0, 1), xticks=[], yticks=[], aspect='equal')
arrows = ax.quiver(*frames[0], angles='xy', scale_units='xy', scale=0.25, width=0.005, pivot='mid')

def draw(i):
    x, y, vx, vy = frames[i]
    arrows.set_offsets(np.c_[x, y])
    arrows.set_UVC(vx, vy)
    ax.set_title(f'Step {3*i}')

anim = animation.FuncAnimation(fig, draw, frames=len(frames), interval=100)
plt.close(fig)
HTML(anim.to_jshtml(default_mode='loop'))
Sim: boids, as 200 agents, steps 0 to 200
  ──────────────────────────────   0%  step 0  0.0 s
  ──────────────────────────────   1%  step 1  0.0 s
  ━─────────────────────────────   2%  step 4  0.0 s
  ━─────────────────────────────   4%  step 7  0.0 s
  ━━────────────────────────────   5%  step 10  0.0 s
  ━━────────────────────────────   7%  step 13  0.0 s
  ━━━───────────────────────────   8%  step 16  0.0 s
  ━━━───────────────────────────  10%  step 19  0.0 s
  ━━━───────────────────────────  11%  step 22  0.0 s
  ━━━━──────────────────────────  13%  step 25  0.0 s
  ━━━━──────────────────────────  14%  step 28  0.0 s
  ━━━━━─────────────────────────  16%  step 31  0.0 s
  ━━━━━─────────────────────────  17%  step 34  0.0 s
  ━━━━━━────────────────────────  19%  step 37  0.0 s
  ━━━━━━────────────────────────  20%  step 40  0.0 s
  ━━━━━━━───────────────────────  22%  step 43  0.0 s
  ━━━━━━━───────────────────────  23%  step 46  0.0 s
  ━━━━━━━───────────────────────  25%  step 49  0.0 s
  ━━━━━━━━──────────────────────  26%  step 52  0.0 s
  ━━━━━━━━──────────────────────  28%  step 55  0.0 s
  ━━━━━━━━━─────────────────────  29%  step 58  0.0 s
  ━━━━━━━━━─────────────────────  31%  step 61  0.0 s
  ━━━━━━━━━━────────────────────  32%  step 64  0.0 s
  ━━━━━━━━━━────────────────────  34%  step 67  0.0 s
  ━━━━━━━━━━━───────────────────  35%  step 70  0.0 s
  ━━━━━━━━━━━───────────────────  37%  step 73  0.0 s
  ━━━━━━━━━━━───────────────────  38%  step 76  0.0 s
  ━━━━━━━━━━━━──────────────────  40%  step 79  0.0 s
  ━━━━━━━━━━━━──────────────────  41%  step 82  0.0 s
  ━━━━━━━━━━━━━─────────────────  43%  step 85  0.0 s
  ━━━━━━━━━━━━━─────────────────  44%  step 88  0.0 s
  ━━━━━━━━━━━━━━────────────────  46%  step 91  0.0 s
  ━━━━━━━━━━━━━━────────────────  47%  step 94  0.0 s
  ━━━━━━━━━━━━━━━───────────────  49%  step 97  0.0 s
  ━━━━━━━━━━━━━━━───────────────  50%  step 100  0.0 s
  ━━━━━━━━━━━━━━━━──────────────  52%  step 103  0.0 s
  ━━━━━━━━━━━━━━━━──────────────  53%  step 106  0.0 s
  ━━━━━━━━━━━━━━━━──────────────  55%  step 109  0.0 s
  ━━━━━━━━━━━━━━━━━─────────────  56%  step 112  0.0 s
  ━━━━━━━━━━━━━━━━━─────────────  58%  step 115  0.0 s
  ━━━━━━━━━━━━━━━━━━────────────  59%  step 118  0.0 s
  ━━━━━━━━━━━━━━━━━━────────────  61%  step 121  0.0 s
  ━━━━━━━━━━━━━━━━━━━───────────  62%  step 124  0.0 s
  ━━━━━━━━━━━━━━━━━━━───────────  64%  step 127  0.0 s
  ━━━━━━━━━━━━━━━━━━━━──────────  65%  step 130  0.0 s
  ━━━━━━━━━━━━━━━━━━━━──────────  67%  step 133  0.0 s
  ━━━━━━━━━━━━━━━━━━━━──────────  68%  step 136  0.0 s
  ━━━━━━━━━━━━━━━━━━━━━─────────  70%  step 139  0.0 s
  ━━━━━━━━━━━━━━━━━━━━━─────────  71%  step 142  0.0 s
  ━━━━━━━━━━━━━━━━━━━━━━────────  73%  step 145  0.0 s
  ━━━━━━━━━━━━━━━━━━━━━━────────  74%  step 148  0.0 s
  ━━━━━━━━━━━━━━━━━━━━━━━───────  76%  step 151  0.0 s
  ━━━━━━━━━━━━━━━━━━━━━━━───────  77%  step 154  0.0 s
  ━━━━━━━━━━━━━━━━━━━━━━━━──────  79%  step 157  0.0 s
  ━━━━━━━━━━━━━━━━━━━━━━━━──────  80%  step 160  0.0 s
  ━━━━━━━━━━━━━━━━━━━━━━━━──────  82%  step 163  0.0 s
  ━━━━━━━━━━━━━━━━━━━━━━━━━─────  83%  step 166  0.0 s
  ━━━━━━━━━━━━━━━━━━━━━━━━━─────  85%  step 169  0.0 s
  ━━━━━━━━━━━━━━━━━━━━━━━━━━────  86%  step 172  0.0 s
  ━━━━━━━━━━━━━━━━━━━━━━━━━━────  88%  step 175  0.0 s
  ━━━━━━━━━━━━━━━━━━━━━━━━━━━───  89%  step 178  0.0 s
  ━━━━━━━━━━━━━━━━━━━━━━━━━━━───  91%  step 181  0.0 s
  ━━━━━━━━━━━━━━━━━━━━━━━━━━━━──  92%  step 184  0.0 s
  ━━━━━━━━━━━━━━━━━━━━━━━━━━━━──  94%  step 187  0.0 s
  ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━─  95%  step 190  0.0 s
  ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━─  97%  step 193  0.0 s
  ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━─  98%  step 196  0.0 s

From random headings at the start, groups form within a few dozen steps, and by the end the boids fly as a few large flocks (which wrap around the edges of the square), with space between their members.

As in the wealth model, each rule is applied to every boid at once, as operations on arrays, rather than by looping over the boids. Working out the distance between every pair of boids takes time and memory in proportion to the square of their number, which is fine for hundreds or a few thousand; for many more, a spatial index (such as SciPy’s cKDTree) would find each boid’s neighbors faster.

Going further