The Boids!¶

Now that we have covered the basics of both NumPy and Matplotlib, we will go through an extended example of using these libaries to run and visualise a simulation of flocking dynamics.

Flocking¶

The aggregate motion of a flock of birds, a herd of land animals, or a school of fish is a beautiful and familiar part of the natural world... The aggregate motion of the simulated flock is created by a distributed behavioral model much like that at work in a natural flock; the birds choose their own course. Each simulated bird is implemented as an independent actor that navigates according to its local perception of the dynamic environment, the laws of simulated physics that rule its motion, and a set of behaviors programmed into it... The aggregate motion of the simulated flock is the result of the dense interaction of the relatively simple behaviors of the individual simulated birds.

— Craig W. Reynolds, "Flocks, Herds, and Schools: A Distributed Behavioral Model", Computer Graphics 21 4 1987, pp 25-34 (see the original paper)

A basic model of flocking behaviour can be derived from three basic rules

• Collision avoidance: avoid collisions with nearby flockmates.
• Velocity matching: attempt to match velocity with nearby flockmates.
• Flock centering: attempt to stay close to nearby flockmates.

We will use a simple two-dimensional model in which the state of the n^th boid (short for bird-oid or flock member) is described by a position vector $\boldsymbol{x}_n = (x_{n,0}, x_{n,1})$ and velocity $\boldsymbol{v}_n = (v_{n,0}, v_{n,1})$. We will represent our flock state as NumPy arrays, implement our simulation dynamics using NumPy array operations and use the animation capabilities of Matplotlib to create animated simulations of our flying boids. We first import the relevant modules numpy, matlotlib.pyplot and matplotlib.animation and set Matplotlib animations to default to being represented as interactive JavaScript widgets in the notebook interface.

Initialising the simulation¶

We will represent the positions and velocities of all the boids as a pair of two-dimensional arrays, both with shape (num_boids, 2) where num_boids is an integer determining the number of boids. We need to set initial values for these arrays, representing the positions and velocities at the start of the simulation. Here we will use a NumPy random number generator object to sample random initial values for the positions and velocities from some distribution.

Generation of random number is a very common task in research software, but a topic with lots of important subtleties. A common pattern in many programming languages is to use a global random number generator which has a state which is updated by every call to random number generation routines. Use of global state like this breaks the usual expectation that for fixed inputs a function will produce fixed outputs, can lead to hard to diagnose bugs and make it difficult to ensure reproducibility of our research.

To avoid these issues we will create an instance of the NumPy Generator class, which encapsulates the state of a random number generator plus methods to produce random numbers from specified distributions, and explicitly pass this in to functions in which we wish to generate random numbers. We will also fix the initial state of the random number generator by explicitly specifying a seed value, meaning we can reproduce our results. The easiest way to produce a seeded Generator object is using the np.random.default_rng function and passing in integer seed argument

The rng object we just created provides methods for producing from various different probability distributions. Here we generate random values uniformly distributed over a specified interval using the Generator.uniform method. We can inspect the docstring for this method by running rng.uniform? (if in a Jupyter notebook or IPython interpreter) or help(rng.uniform) (in any Python interpreter).

Helpfully, Generator.uniform accepts a size argument allowing us to generate an array of random samples of a specified shape (the argument is named size rather than shape to avoid clashing with the usage of shape as the standard name for a parameter of some probability distributions). We can also specify array-like arguments for the low and high parameters specifying the lower and upper bounds of the interval the distribution is defined on, these should either match the shape specified by size or be of a compatible shape for broadcasting. Below we define a function which given a Generator object rng and optional arguments specifying the number of boids num_boids, minimum and maximum of the intervals to draw the positions (min_position and max_position) and velocities (min_velocity and max_velocity) from, outputs a pair of arrays positions and velocities, both of shape (num_boids, 2), corresponding to the sampled position and velocities values respectively.

We can call this function with our rng random number generator object (and default values for the other arguments) to generate random initial values for the positions and velocities arrays.

Visualising the boids¶

Now that we have initialised the state of our simulation, we are ready to start visualising our boids. As we have assumed our boids exist in two-dimensions, we can use Matplolib's extensive range of plotting functions for two-dimensional data. Ideally we want to represent both the instantenous positions and velocities of the boids in our visualisation. We will use a Matplotlib quiver plot to do this, representing each boid as an arrow located on the plot axes at a point corresponding to its simulated position, and with the direction of the arrow representing its current heading. While we could also represent its speed by the size of the arrow, we choose here to keep the arrows all the same size - this requires calculating the unit vectors (vectors of length one) corresponding to the velocity vectors, which can be done efficiently with broadcasted array operations in NumPy. We wrap the code for plotting the boids into a function to allow for easy reuse. We return the generated Matplotlib figure and axis (subplot) objects, along with artist object representing the quiver plot arrows returned by the quiver call, to allow us to make calls to these objects outside of the function.

We can now produce a (rather boring) visualisation of our boids in their initial state.

Simulating the model dynamics¶

We will model the boids as being subject to forces corresponding to the three rules described above. Newton's second law of motion tells use that the acceleration of a body (rate of change of velocity) is proportional to the net force acting on it, with the velocity in turn being the rate in change of the position of the object. In mathematical notation, if $F_{k,n}$ is a function which calculates the value of the $k$th force on the $n$th boid given the positions and velocities of all the boids then we have that

$$ \frac{\mathrm{d}\boldsymbol{x}_{0:N}}{\mathrm{d}t} = \boldsymbol{v}_{0:N}, \qquad \frac{\mathrm{d}\boldsymbol{v}_{0:N}}{\mathrm{d}t} = \sum_{k=1}^K F_{k,n}(\boldsymbol{x}_{0:N}, \boldsymbol{v}_{0:N}). $$

Note that it is not a problem if you are not familiar with the notation here - understanding the mathematics behind the model is non-essential! There are a variety of numerical approaches for approximately simulating dynamics of this form. Here we will use a particularly simple to implement numerical method which corresonds to the approximation that for a small timestep $\delta t$ that

$$ \boldsymbol{x}_{0:N}(t + \delta t) \approx \boldsymbol{x}_{0:N}(t) + \delta t \boldsymbol{v}_{0:N}(t), \qquad \boldsymbol{v}_{0:N}(t + \delta t) \approx \boldsymbol{x}_{0:N}(t) + \delta t \sum_{k=1}^K F_k(\boldsymbol{x}_{0:N}(t + \delta t), \boldsymbol{v}_{0:N}(t)). $$

Using NumPy array operations we can implement this numerical scheme efficiently by updating the positions and velocities for all boids simulateneously. We wrap this into a function simulate_timestep below that updates the positions and velocities arrays in-place.

Creating an animation of the simulation¶

Now that we have a function to update the boid states according to the model dynamics, we are ready to produce animations visualising the simulated dynamics over time. To do this we will use the FuncAnimation class from the matplotlib.animation module. We can inspect the docstring for the initialiser for this class:

We see that FuncAnimation produces an animation by repeatedly calling a specified function to update the animation 'frame', which corresponds to a Matplotlib figure instance. As well as the figure object to use and function updating the frames, we need to specify the frames argument - here we do this with an integer corresponding to the number of frames to animate. Finally we specify an optional argument interval corresponding to the delay between each frame being shown in milliseconds. Again we wrap all of this into a function to allow for easy reuse.

We are now ready to produce our first animation! We generate initial positions and velocities using our initialise_boid_states function and then pass these to the animate_flock function, using the default values for the other arguments for now, with in particular we not specifying any forces acting on the boids for now.

Though our boids are now moving, we see that they have the rather uninteresting behaviour of moving perpetually in a straight line. To get more interesting flocking like behaviour we need to specify some forces.

Specifying the flocking dynamics¶

Now that we have our simulation and animation framework set up, we are ready to start defining the forces corresponding to the flocking behaviour rules we encountered at the start of the notebook.

Cohesion: staying close to nearby flockmates¶

Credit: Craig Reynolds (public domain)

The first behaviour we consider is the tendency for flocks to remain close to flockmates, which we term cohesion. We represent this as a force which pushes the flock members towards the mean point of surrounding flock members, that is the force exerted is proportional to the difference between the mean flock position and the current flock members position. We can implement this efficiently using broadcasted NumPy array operations:

Now that we have defined our first force function, we can create a new animation with this force applied:

Our boids now show a more interesting dynamic behaviour, staying together as one cohesive unit. We see however that a lot of the time boids appear to be in (near) collisions with each other.

Separation: avoiding collisions with nearby flockmates¶

Credit: Craig Reynolds (public domain)

The second behaviour that we consider it tendency of flockmates to avoid collisions with each other. We represent this as a force which for each pair of boids within a certain distance of each other, exerts a force which pushes the boids away from each other so that the displacement beween them increases, by acting along the line of displacement. To only sum the displacements over the pairs of boids within a specified distance of each other we can use the numpy.where function. Calling

values = np.where(conditions, values_if_true, values_if_false)

creates an array values such that for an index i (either a single integer for one-dimensional arrays or tuple of integers for multidimensional arrays) values[i] is equal to values_if_true[i] if conditions[i] is True and values_if_false[i] otherwise. All of the array arguments conditions, values_if_true and values_if_false should be of compatible shapes - that is of the same shape, or of shapes that can be broadcast (remember that a scalar is compatible for broadcasting with an array of any shape).

We can now run and animate a simulation with both the forces we have defined:

Alignment: matching velocity with nearby flockmates¶

Credit: Craig Reynolds (public domain)

The final rule we consider is a little different in that it specifies a relationship between the velocities rather than positions of the boids. The alignment rule specifies that nearby flockmates should tend to match velocities with each other. We can implement this as a force which is negatively proportional to the velocity differences between pairs of boids within a certain distance of each other. We can use a similar implementation to that for separation_force:

Now we are finally ready to run and animate a simulation with forces for all three rules defined:

Conclusion and extensions¶

Hopefully this example has illustrated the power and convenience of NumPy. Implementing an equivalent model using native Python data structures such as lists to represent the states of the boids would both result in code that was slower to run and harder to read.

If you are interested in exploring this model further here are some ideas for extensions

• Currently the cohesion force applies globally rather than only acting locally on boids within a certain distance of each other like the other two forces. Can you alter the implementation to have a similar behaviour of only acting over a finite distance?
• There are various model parameters such as the relative force strengths cohesion_strength, separation_strength and alignment_strength which are currently left as the default values specified in the force function definitions. Can you redesign the simulate_timestep and/or animate_flock to allow passing in a dictionary of these parameter values to override the default to allow easily visualising the behaviour for different parameters?
• Currently our boids exist in a two-dimensional world - can you alter the simulation to work in three-dimensions? You may find Matplotlib's three-dimensional plotting toolkit useful for the visualisation side.