If you have ever looked up at the sky in early evening you might have been greeted by the spectacle of many birds flying in unison creating elaborate and complex patterns in the sky. Or perhaps on the discovery channel you have seen a documentary on the sea and have seen schools of fish creating similar elaborate patterns.

There are many evolutionary reasons for why these patterns may form: firstly is the "safety in numbers" concept whereby having many individuals grouped together protects against predators. Predators can only attack the edge of the pattern, the number of individuals on the edges is significantly fewer than the individuals in the middle of the pattern - thus increasing safety for the largest number of individuals. The patterns may also trick predators into thinking the pattern is some larger threatening entity and so less likely to attack in the first place. There is also the possiblity that by creating a group pattern that is large enough to be seen from far away it attracts more individuals to join the group.

This phenomena is known as "flocking" - in this blog post we will investigate how flocks form.

## An Agent Based Approach

One way in which we could reverse engineer these patterns is to have a "central controller" that tells where each flock member should be located at each moment in time. While this may produce results that are seemingly indistinguishable from real-world flocking behaviour we would be very hard pressed to argue that this is what happens in nature. Instead we would like each member to have some level of autonomy and a simple limited set of rules that when combined produces the desired result - if such a model can be found it will allow us to analyse and study flock formation in more detail. Any model containing many members each following simple rules within an environment can be classified as an "agent based model" (members herein being referred to as agents). Agent based models span further than flocking models, they can be used anywhere where it's possible to define simple rules for local behaviour and where large scale/population models (e.g. differential equations) do not provide details on the scales that are of interest to us. Agent based models have been used to model a whole manner of areas from biologically inspired predator/prey interactions, spreading of infection through to the social sciences and idea spreading and the economy.

## Model Design

Flocking models have been studied fairly extensively. In this blog I will present my implementation of a flocking model, there will be "prey" agents who are large in number and form flocks to try and protect themselves. There will be a small number of "predator" agents who try and eat/kill the prey. Each class of agent has their own dynamic rules for how they move in space.

### Rules for Prey

For the prey agents we will assume that there are 5 forces operating to direct the agents. These forces are equivalent to acceleration vectors, we will integrate over time to change the velocity and integrate velocity to get the change in position. We can describe the forces at bay as follows:

• Flocking Forces - This represents the attraction of an agent to move towards a density of other agents (join the flock). To do this we consider the velocity vector of agents which we denote $\mathbf{v_i}$. We assume that we will only flock towards others if they are near us, we thus have a "flocking radius" which we denote $r_f$ which is how far the agent can "see". We define a vector quantity: $\mathbf{V_i} = \sum_{j \in N_{r_f}(i)} \mathbf{v_j}$. Where we denote the set of agents within a neighbourhood of radius $r$ of agent $i$ as: $N_{r}(i)$. The strength of this force is defined by a parameter $\alpha$. We then define the flocking force applied to agent $i$ as: $$\mathbf{F}_{flock}^i = \alpha \frac{\mathbf{V_i}}{|\mathbf{V_i}|}$$
• Repulsive Forces - If we are in a crowd we like our own space, if somebody gets to close we move away. At the most extreme level there is a limit on how close one can get to another owing to the size of their bodies. The repulsive force acts to stop agents getting too close together. This force operates only on a radius $r_0$ - if somebody is very far away they won't impact how we move in space. We denote $|\mathbf{d_{ij}}|$ to be the scalar distance (magnitude) between agents $i$ and $j$. We use boldface $\mathbf{d_{ij}}$ to represent the vector distance. The strength of this force is defined by the parameter $\epsilon$. We can then express this force applied to agent $i$ as: $$\mathbf{F}_{rep}^i = \epsilon \sum_{j \in N_{r_0}(i)}\left(1 - \frac{|\mathbf{d_{ij}}|}{r_0}\right)^{\frac{3}{2}} \frac{\mathbf{d_{ij}}}{|\mathbf{d_{ij}}|}$$
• Self-Propulsion Forces - This force represents a desire to "keep moving". In some situations there is a tendency for sustained movement - for example in a march there is a desire to keep moving at a constant velocity. We represent this by the self-propulsion force. The strength of this force is controlled via a parameter $\mu$. For agent $i$ we can express this as: $$\mathbf{F}_{prop}^i = \mu \left(v_0 - |\mathbf{v_i}|\right)\frac{\mathbf{v_i}}{|\mathbf{v_i}|}$$
• Random Forces - In the model we may wish to allow for random forces, that is each agent has it's own autonomy. To do this we will just apply a uniform random draw. We will use a parameter $r_{amp}$ (random amplitude) to denote the strength of this force.
• Predator-Repulsive Forces - One of the driving forces for any prey is to escape predators. In this model we will assume that this force trumps all others, that is if a predator is close an agent will run away and ignore all other desires of flocking/self-propulsion/etc. This operates in the exact same way as the repulsive force above but instead of repulsion from other agents it considers repulsion from predators. We define the radius over which this occurs as $r_p$ and the strength of this force $\delta$ - these replace $r_0$ and $\epsilon$ in the formula above. We also apply a "speed limit" to the agents to prevent them travelling arbitrarily quickly.

### Rules for Predators

For predators the dynamics are much simpler. We assume that a predator has infinite vision, it will look for the nearest prey and run towards it at a fixed speed. However we will add a (possibly small) random component to the movement - this could represent indecisiveness on the part of the predator or the evasive tactics employed by the prey. The only model assumptions for predators is the predator speed $pv_{max}$ and the proportion of movement that is random $pr_{amp}$.

### Environment

We will assume a simple environment consisting of a unit-square. The boundaries will be periodic (i.e. the left edge will "connect" to the right edge and top to bottom). Like the game asteroids this assumes the dynamics take place on the surface of a torus. This assumption simplifies the dynamics since there is no need to consider any boundary conditions - if we assume a fixed "wall" around the environment we would need to define rules as to what happens when an agent tries to cross (e.g. do they "reflect" off walls? Do they "stick"? Do they slow down and change direction to avoid hitting the wall all together?). For convenience in calculating distances in the periodic boundary world we have a helper function called "buffer" which replicates any agents near the walls and places them outside the unit-square - this ensures that distances are calculated correctly.

### Other Model Assumptions

There are a few other modelling assumptions made, one of the major decisions of any agent based model is deciding whether to use synchronous or asynchronous updating. We will assume synchronous updating here, each agent will make their movements at the same time. This avoids the problem of agents who are updated "later" within a tick of an asynchronous scheme gaining an "advantage" over those updated first (this would be a particular issue for predator-prey relations where if a predator waits until the prey has moved it is more likely to catch the prey!)

To avoid the problem of a predator following a single prey-agent aimlessly we will use an "eat" dynamic. Here if a predator gets within a small distance of a prey it will "eat" it, the prey will disappear and instantly re-spawn to a random location in the environment with a random starting velocity. There will be no limit to how many prey a predator can eat in a single step and it is assumed the predator can eat without stopping/slowing down.

## The Code

An implementation of this model can be seen below. I use matplotlib's animaton feature to make an animation of the system through time. In the animation prey are represented by small blue crosses and predators by larger red diamonds.

import numpy as np
from matplotlib import animation
import matplotlib.pyplot as plt
from IPython.display import HTML
%matplotlib inline

###### Model Parameters ######
dt = 0.01           # Time step
num_frames = 250    # Number of animated frames
burn = 1000         # Number of burn in iterations

# Flock Parameters
N = 250             # Number of flocking agents
r0 = 0.1            # Repulsion force range
eps = 2.0           # Repulsion force strength
rf = 0.25           # Flocking force range
alpha = 10.0         # Flocking force strength
v0 = 0.25           # Target speed
vmax = 0.75         # Maximum speed
mu = 0.1            # Self-Propulsion force
ramp = 1.0          # Random force
rp = 0.5            # Predator-repulsion force range
delta = 5.0         # Predator-repulsion force strength

# Predator Parameters
M = 3               # Number of predator agents
pv_max = 1.0        # Predator speed
pramp = 0.25         # Predator random force
eat = 0.025         # Eating radius

###### Fix Random Seed ######
seed = 123
np.random.seed(seed)

###### Helper Functions ######
def buffer(rb, x, y, vx, vy):
"""
buffer(rb, x, y, vx, vy)

Takes copies of agents near the edges
of the unit square and extends outwards
by a buffer margin rb. This makes it
easier to calculate distances between
agents for flocking and repulsive forces

Inputs:
rb - buffer range
x - array of x coordinate positions [0,1]
y - array of y coordinate positions [0,1]
vx - array of x-component of velocity
vy - array of y-component of velocity

Outputs:
nb, xb, yb, vxb, vyb - Buffered copies of inputs
"""
_N = x.shape[0]

# Initialize buffer arrays
xb[0:_N] = x
yb[0:_N] = y
vxb[0:_N] = vx
vyb[0:_N] = vy

# nb is a counter already have _N agents
nb = _N-1

for i in range(nb+1):
if x[i] <= rb:
# left edge
nb+=1
xb[nb] = x[i]+1
yb[nb] = y[i]
vxb[nb] = vx[i]
vyb[nb] = vy[i]
if x[i] >= 1 - rb:
# right edge
nb+=1
xb[nb] = x[i]-1
yb[nb] = y[i]
vxb[nb] = vx[i]
vyb[nb] = vy[i]
if y[i] <= rb:
# bottom edge
nb+=1
xb[nb] = x[i]
yb[nb] = y[i]+1
vxb[nb] = vx[i]
vyb[nb] = vy[i]
if y[i] >= 1 - rb:
# top edge
nb+=1
xb[nb] = x[i]
yb[nb] = y[i]-1
vxb[nb] = vx[i]
vyb[nb] = vy[i]
if (x[i] <= rb) and (y[i] <= rb):
# bottom left corner
nb+=1
xb[nb] = x[i]+1
yb[nb] = y[i]+1
vxb[nb] = vx[i]
vyb[nb] = vy[i]
if (x[i] <= rb) and (y[i] >= 1 - rb):
# top left corner
nb+=1
xb[nb] = x[i]+1
yb[nb] = y[i]-1
vxb[nb] = vx[i]
vyb[nb] = vy[i]
if (x[i] >= 1 - rb) and (y[i] <= rb):
# bottom right corner
nb+=1
xb[nb] = x[i]-1
yb[nb] = y[i]+1
vxb[nb] = vx[i]
vyb[nb] = vy[i]
if (x[i] >= 1 - rb) and (y[i] >= 1 - rb):
# top right corner
nb+=1
xb[nb] = x[i]-1
yb[nb] = y[i]-1
vxb[nb] = vx[i]
vyb[nb] = vy[i]

return nb, xb, yb, vxb, vyb

def force(nb, xb, yb, vxb, vyb, x, y, vx, vy, mb, pxb, pyb):
"""
force(nb, xb, yb, vxb, vyb, x, y, vx, vy, mb, pxb, pyb)

Calculate force applied to agents (x,y) divided into:
- flocking force (flockx, flocky)
- repulsive force (repx, repy)
- Self propulsion (fpropx, fpropy)
- Random force (frandx, frandy)
- Predator repulsive force (fpredx, fpredy)

Inputs:
nb - number of buffered agents
xb - buffered agent x position array
yb - buffered agent y position array
vxb - buffered agent velocity x coordinate array
vyb - buffered agent velocity y coordinate array
x - agent x position array
y - agent y position array
vx - agent x coordinate velocity array
vy - agent y coordinate velocity array
mb - number buffered predator agents
pxb - predator x-coordinate position array
pyb - predator y-coordinate position array

Global Variables Called:
rf - Flocking force range
r0 - Repulsion force range
alpha - Flocking force strength
mu - Self propulsion force
v0 - Target velocity
ramp - Random force strength
rp - Predator repulsion force range

Outputs:
fx, fy - x and y coordinate forces
"""

_N = x.shape[0]

for i in range(_N):
repx = 0
repy = 0
flockx = 0
flocky = 0
nflock = 0
fpredx = 0
fpredy = 0

for j in range(nb):
d2 =  (xb[j] - x[i])**2 + (yb[j] - y[i])**2

if (d2 <= rf**2) and (i != j):
flockx += vxb[j]
flocky += vyb[j]
nflock += 1

if d2 <= r0**2:
d = np.sqrt(d2)
repx += eps*(1- d / r0)**1.5 * (x[i] - xb[j]) / d
repy += eps*(1- d / r0)**1.5 * (y[i] - yb[j]) / d

normflock = np.sqrt(flockx**2 + flocky**2)
if nflock == 0:
normflock = 1
flockx = alpha * flockx / normflock
flocky = alpha * flocky / normflock
vnorm = np.sqrt(vx[i]**2 + vy[i]**2)
fpropx = mu * (v0 - vnorm) * (vx[i] / vnorm)
fpropy = mu * (v0 - vnorm) * (vy[i] / vnorm)
frandx = ramp * (2*np.random.random() - 1)
frandy = ramp * (2*np.random.random() - 1)
fx[i] = (flockx + frandx + fpropx + repx)
fy[i] = (flocky + frandy + fpropy + repy)

for k in range(mb):
d2 =  (pxb[k] - x[i])**2 + (pyb[k] - y[i])**2

if (d2 <= rp**2):
d = np.sqrt(d2)
fpredx += delta*(1- d / rp)**1.5 * (x[i] - pxb[k]) / d
fpredy += delta*(1- d / rp)**1.5 * (y[i] - pyb[k]) / d

if fpredx**2 + fpredy**2 > 0:
fx[i] = fpredx
fy[i] = fpredy

return fx, fy

def pred_dist(px1, py1, x1, y1):
"""
pred_dist(px1, py1, x1, y1)

Returns an array of distances from predator to an array of prey

Inputs:
px1 - predator x coordinate
py1 - predator y coordinate
x1 - array of prey x coordinates
y1 - array of prey y coordinates

Outputs:
distance array
"""
dx1 = np.abs(px1-x1)
dx2 = np.abs(1+px1-x1)
dx3 = np.abs(px1-(1+x1))
dy1 = np.abs(py1-y1)
dy2 = np.abs(1+py1-y1)
dy3 = np.abs(py1-(1+y1))

dx = np.minimum(np.minimum(dx1, dx2), dx3)
dy = np.minimum(np.minimum(dy1, dy2), dy3)

return np.sqrt(dx**2+dy**2)

def direction(start_x, start_y, target_x, target_y):
"""
direction(start_x, start_y, target_x, target_y)

Find a unit direction vector from a starting point to a target

Inputs:
start_x - starting point x coordinate
start_y - starting point y coordinate
target_x - target point x coordinate
target_y - target point y coordinate

Returns
dir_x - direction in x
dir_y - direction in y
"""
dx1 = target_x - start_x
dx2 = target_x - (start_x + 1)
dx3 = target_x - (start_x - 1)
dy1 = target_y - start_y
dy2 = target_y - (start_y + 1)
dy3 = target_y - (start_y - 1)

if abs(dx1) <= (abs(dx2) and abs(dx3)):
dir_x = dx1
elif abs(dx2) <= (abs(dx1) and abs(dx3)):
dir_x = dx2
elif abs(dx3) <= (abs(dx1) and abs(dx2)):
dir_x = dx3
if abs(dy1) <= (abs(dy2) and abs(dy3)):
dir_y = dy1
elif abs(dy2) <= (abs(dy1) and abs(dy3)):
dir_y = dy2
elif abs(dy3) <= (abs(dy1) and abs(dy2)):
dir_y = dy3

return dir_x, dir_y

def unitize(vec_x, vec_y):
"""
unitize(vec_x, vec_y)

Returns a vector of same direction as input
but with unit size
"""
norm = np.sqrt(vec_x**2 + vec_y**2)

return vec_x/norm, vec_y/norm

def update():
"""
update()

Update function that updates position and
velocities of prey and predator agents by
one time step

Inputs:
None

Outputs:
None (updates values of global variables)
"""

global x
global y
global vx
global vy
global px
global py
global vpx
global vpy
global eat_count

x_old = x.copy()
y_old = y.copy()

nb, xb, yb, vxb, vyb = buffer(max(r0,rf), x, y, vx, vy)
mb, pxb, pyb, vpxb, vpyb = buffer(rp, px, py, vpx, vpy)
fx, fy = force(nb, xb, xb, vxb, vyb, x, y, vx, vy, mb, pxb, pyb)

# Use force to calculate speeds
vx += fx * dt
vy += fy * dt

# Apply speed limit of maximum velocity
vmod = np.sqrt(vx**2 + vy**2)
vmult = np.where(vmod < vmax, vmod, vmax) / vmod

vx *= vmult
vy *= vmult

# Calculate new positions and re-position in unit square
x += vx * dt
y += vy * dt
x = (1 + x) % 1
y = (1 + y) % 1

# Predator calculation
for i in range(M):
pred_x = px[i]
pred_y = py[i]

# Calculate distances to prey
dist = pred_dist(pred_x, pred_y, x_old, y_old)

# If prey within eat range then eat and regenerate prey
for j in range(N):
if dist[j] < eat:
x[j] = np.random.random()
y[j] = np.random.random()
vx[j] = 2*np.random.random() - 1
vy[j] = 2*np.random.random() - 1
eat_count += 1

# Find closest prey
min_dist = np.min(dist)

# Calculate direction to move in
vpx_t, vpy_t = direction(pred_x, pred_y, prey_x, prey_y)
vpx_t, vpy_t = unitize(vpx_t, vpy_t)

# Simulate random component and normalize
vpx_rand = 2*np.random.random() - 1
vpy_rand = 2*np.random.random() - 1
vpx_rand, vpy_rand = unitize(vpx_rand, vpy_rand)

# Combine prey-direction and random components
vpx[i] = (1-pramp)*vpx_t + pramp*vpx_rand
vpy[i] = (1-pramp)*vpy_t + pramp*vpy_rand
vpx[i], vpy[i] = unitize(vpx[i], vpy[i])

# Scale to maximum velocity
vpx[i] *= pv_max
vpy[i] *= pv_max

# Move predator forward
px[i] += vpx[i]*dt
py[i] += vpy[i]*dt

# Re-position predators to unit square
px = (1 + px) % 1
py = (1 + py) % 1

###### Main Code ######

###### Set Up Arrays ######
# Flock arrays
x = np.zeros(N)
y = np.zeros(N)
vx = np.zeros(N)
vy = np.zeros(N)
fx = np.zeros(N)
fy = np.zeros(N)
# Flock buffer arrays
xb = np.zeros([4*N])
yb = np.zeros([4*N])
vxb = np.zeros([4*N])
vyb = np.zeros([4*N])

# Predator arrays
px = np.zeros(M)
py = np.zeros(M)
vpx = np.zeros(M)
vpy = np.zeros(M)
# Predator buffer arrays
pxb = np.zeros([4*M])
pyb = np.zeros([4*M])
vpxb = np.zeros([4*M])
vpyb = np.zeros([4*M])

# Eat counter initialization
eat_count = 0

# Initialize positions and velocities
for i in range(N):
x[i] = np.random.random()
y[i] = np.random.random()
vx[i] = 2*np.random.random() - 1
vy[i] = 2*np.random.random() - 1

for i in range(M):
px[i] = np.random.random()
py[i] = np.random.random()
vpx[i] = 2*np.random.random() - 1
vpy[i] = 2*np.random.random() - 1

# Burn simulations to improve initial conditions
for i in range(burn):
update()

# Reset Eat count
eat_count = 0

###### Animate ######
# Initialize figure
figure = plt.figure()
axes = plt.axes(xlim=(0, 1), ylim=(0, 1))
axes.set_xticks([], [])
axes.set_yticks([], [])
scatter_prey = axes.scatter(x, y, marker='x', s=25, c='Blue')
scatter_pred = axes.scatter(px, py, marker='D', s=100, c='Red')

# Define animation step
def animate(frame):
update()
prey_data =  np.array((x, y)).T
scatter_prey.set_offsets(prey_data)
pred_data =  np.array((px, py)).T
scatter_pred.set_offsets(pred_data)

# Display animation function
def display_animation(anim):
# Close initialized static plot
plt.close(anim._fig)

# Returns animation to HTML
return HTML(anim.to_jshtml())

# Set up FuncAnimation
anim = animation.FuncAnimation(figure, animate, frames=num_frames, interval=50)

# Call the display function
display_animation(anim)

</input>