{
"cells": [
{
"cell_type": "markdown",
"id": "bc8120df",
"metadata": {},
"source": [
"# The Boids!\n",
"\n",
"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.\n",
"\n",
"## Flocking\n",
"\n",
"\n",
"> The aggregate motion of a flock of birds, a herd of land animals, or a school of fish is a beautiful and familiar\n",
"part of the natural world... The aggregate motion of the simulated flock is created by a distributed behavioral model much\n",
"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. \n",
"\n",
"— Craig W. Reynolds, \"Flocks, Herds, and Schools: A Distributed Behavioral Model\", *Computer Graphics* **21** _4_ 1987, pp 25-34 (see the [original paper](http://www.cs.toronto.edu/~dt/siggraph97-course/cwr87/))\n",
"\n",
"A basic model of flocking behaviour can be derived from [three basic rules](https://en.wikipedia.org/wiki/Boids)\n",
"\n",
"> * *Collision avoidance*: avoid collisions with nearby flockmates.\n",
"> * *Velocity matching*: attempt to match velocity with nearby flockmates.\n",
"> * *Flock centering*: attempt to stay close to nearby flockmates.\n",
"\n",
"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](https://matplotlib.org/stable/api/animation_api.html) 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](http://louistiao.me/posts/notebooks/embedding-matplotlib-animations-in-jupyter-as-interactive-javascript-widgets/)."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "e5a8cbf6",
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"from matplotlib import pyplot as plt, animation\n",
"plt.rcParams['animation.html'] = 'jshtml'"
]
},
{
"cell_type": "markdown",
"id": "7c146906",
"metadata": {},
"source": [
"## Initialising the simulation\n",
"\n",
"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. \n",
"\n",
"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.\n",
"\n",
"To avoid these issues we will create an instance of [the NumPy `Generator` class](https://numpy.org/doc/stable/reference/random/generator.html), 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"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "d1cf77f2",
"metadata": {},
"outputs": [],
"source": [
"rng = np.random.default_rng(seed=21878533808081494313)"
]
},
{
"cell_type": "markdown",
"id": "bd65e260",
"metadata": {},
"source": [
"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](https://numpy.org/doc/stable/reference/random/generator.html). 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).\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "b930efbf",
"metadata": {},
"outputs": [],
"source": [
"rng.uniform?"
]
},
{
"cell_type": "markdown",
"id": "2cd5b289",
"metadata": {},
"source": [
"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."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "3f194d4a",
"metadata": {},
"outputs": [],
"source": [
"def initialise_boid_states(\n",
" rng, \n",
" num_boids=100,\n",
" min_position=(100, 900), \n",
" max_position=(200, 1100), \n",
" min_velocity=(0, -20), \n",
" max_velocity=(10, 20)\n",
"):\n",
" \"\"\"Generate random initial states for the boids.\n",
" \n",
" Args:\n",
" rng: NumPy random number generator object.\n",
" num_boids: Number of boids to generate states for.\n",
" min_position: Length 2 sequence defining lower bounds for \n",
" interval to uniformly generate positions from.\n",
" max_position: Length 2 sequence defining upper bounds for\n",
" interval to uniformly generate position from.\n",
" min_velocity: Length 2 sequence defining lower bounds for \n",
" interval to uniformly generate velocities from.\n",
" max_velocity: Length 2 sequence defining upper bounds for\n",
" interval to uniformly generate velocities from.\n",
" \n",
" Returns:\n",
" Tuple of two arrays `(positions, velocities)`, both of shape\n",
" `(num_boids, 2)` corresponding to respectively the positions and\n",
" velocities of boids.\n",
" \"\"\"\n",
" positions = rng.uniform(min_position, max_position, size=(num_boids, 2))\n",
" velocities = rng.uniform(min_velocity, max_velocity, size=(num_boids, 2))\n",
" return positions, velocities"
]
},
{
"cell_type": "markdown",
"id": "5b2f280a",
"metadata": {},
"source": [
"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."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "2cf17a75",
"metadata": {},
"outputs": [],
"source": [
"positions, velocities = initialise_boid_states(rng)"
]
},
{
"cell_type": "markdown",
"id": "ce61f5e0",
"metadata": {},
"source": [
"## Visualising the boids\n",
"\n",
"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](https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.quiver.html) 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."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "cd48009e",
"metadata": {},
"outputs": [],
"source": [
"def calculate_unit_vectors(vectors):\n",
" \"\"\"Calculate the length-one vectors corresponding to an array of vectors.\n",
" \n",
" Args:\n",
" vectors: Array with last dimension corresponding to vector dimension.\n",
" \n",
" Returns:\n",
" Array of same shape as `vectors` with Euclidean norm along last axis\n",
" equal to one.\n",
" \"\"\"\n",
" return vectors / (vectors**2).sum(-1)[..., np.newaxis]**0.5\n",
"\n",
"\n",
"def plot_boids(positions, velocities, figsize=(8, 8), xlim=(0, 2000), ylim=(0, 2000)):\n",
" \"\"\"Create visual representation of boids as quiver plot.\n",
" \n",
" Args:\n",
" positions: Array of shape `(num_boids, 2)` defining positions of boids.\n",
" velocities: Array of shape `(num_boids, 2)` defining velocities of boids.\n",
" figsize: Tuple defining figure dimension `(width, height)` in inches.\n",
" xlim: Tuple `(min, max)` defining extents of horizontal axis.\n",
" ylim: Tuple `(min, max)` defining extents of vertical axis.\n",
" \n",
" Returns:\n",
" Tuple containing Matplotlib figure, axis and artist corresponding to plot.\n",
" \"\"\"\n",
" fig, ax = plt.subplots(figsize=(8, 8))\n",
" ax.set(xlim=xlim, ylim=ylim, xlabel=\"$x_0$ coordinate\", ylabel=\"$x_1$ coordinate\")\n",
" velocity_unit_vectors = calculate_unit_vectors(velocities)\n",
" arrows = ax.quiver(\n",
" positions[:, 0], # horizontal coordinates for origin points for arrows\n",
" positions[:, 1], # vertical coordinates for origin points of arrows\n",
" velocity_unit_vectors[:, 0], # horizontal component of arrow vectors\n",
" velocity_unit_vectors[:, 1], # vertical component of arrow vectors\n",
" scale=40, # size of arrows\n",
" angles='xy', # convention used for specifying arrow directions\n",
" color='C0', # color of arrows - set to first value in default color cycle\n",
" pivot='middle' # plot middle of arrows at specified origin points\n",
" )\n",
" return fig, ax, arrows"
]
},
{
"cell_type": "markdown",
"id": "c2d48f0c",
"metadata": {},
"source": [
"We can now produce a (rather boring) visualisation of our boids in their initial state."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "a01dd8b3",
"metadata": {},
"outputs": [],
"source": [
"fig, ax, arrows = plot_boids(positions, velocities)"
]
},
{
"cell_type": "markdown",
"id": "45046e64",
"metadata": {},
"source": [
"## Simulating the model dynamics\n",
"\n",
"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\n",
"\n",
"\n",
"$$\n",
" \\frac{\\mathrm{d}\\boldsymbol{x}_{0:N}}{\\mathrm{d}t} = \\boldsymbol{v}_{0:N}, \\qquad\n",
" \\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}).\n",
"$$\n",
"\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\n",
"\n",
"$$\n",
" \\boldsymbol{x}_{0:N}(t + \\delta t) \\approx \\boldsymbol{x}_{0:N}(t) + \\delta t \\boldsymbol{v}_{0:N}(t), \\qquad\n",
" \\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)).\n",
"$$\n",
"\n",
"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."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "4023ccef",
"metadata": {},
"outputs": [],
"source": [
"def simulate_timestep(positions, velocities, forces, timestep):\n",
" \"\"\"Simulate model dynamics forward one timestep, updating states in-place.\n",
" \n",
" Args:\n",
" positions: Array of shape `(num_boids, 2)` defining positions of boids.\n",
" velocities: Array of shape `(num_boids, 2)` defining velocities of boids.\n",
" forces: Sequence of functions computing forces on boids given positions\n",
" and velocities of all boids.\n",
" timestep: Scalar timestep to use for numerical integrator.\n",
" \"\"\"\n",
" positions += timestep * velocities\n",
" velocities += timestep * sum(force(positions, velocities) for force in forces)"
]
},
{
"cell_type": "markdown",
"id": "82d55a3a",
"metadata": {},
"source": [
"## Creating an animation of the simulation\n",
"\n",
"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](https://matplotlib.org/stable/api/_as_gen/matplotlib.animation.FuncAnimation.html). We can inspect the docstring for the initialiser for this class:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "9ebc411c",
"metadata": {},
"outputs": [],
"source": [
"animation.FuncAnimation?"
]
},
{
"cell_type": "markdown",
"id": "31070c88",
"metadata": {},
"source": [
"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."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "700b5b4f",
"metadata": {},
"outputs": [],
"source": [
"def animate_flock(positions, velocities, forces=(), timestep=1., num_step=100):\n",
" \"\"\"Visualise the dynamics of the boids as a Matplotlib animation.\n",
" \n",
" Args:\n",
" positions: Array of shape `(num_boids, 2)` defining positions of boids.\n",
" velocities: Array of shape `(num_boids, 2)` defining velocities of boids.\n",
" forces: Sequence of functions computing forces on boids given positions\n",
" and velocities of all boids.\n",
" timestep: Scalar timestep to use for numerical integrator.\n",
" num_step: Number of timesteps to simulate in animation.\n",
" \n",
" Returns:\n",
" Matplotlib animation of simulated boid dynamics.\n",
" \"\"\"\n",
" fig, ax, arrows = plot_boids(positions, velocities)\n",
"\n",
" def update_frame(frame_index):\n",
" simulate_timestep(positions, velocities, forces, timestep)\n",
" velocity_unit_vectors = calculate_unit_vectors(velocities)\n",
" arrows.set_offsets(positions)\n",
" arrows.set_UVC(velocity_unit_vectors[:, 0], velocity_unit_vectors[:, 1])\n",
" return [arrows]\n",
" \n",
" # Close Matplotlib figure object to avoid displaying static figure as well as animation\n",
" plt.close(fig)\n",
" return animation.FuncAnimation(fig, update_frame, num_step, interval=50)"
]
},
{
"cell_type": "markdown",
"id": "b437471c",
"metadata": {},
"source": [
"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."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "3a67d379",
"metadata": {},
"outputs": [],
"source": [
"positions, velocities = initialise_boid_states(rng)\n",
"animate_flock(positions, velocities)"
]
},
{
"cell_type": "markdown",
"id": "63abb7c0",
"metadata": {},
"source": [
"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.\n",
"\n",
"## Specifying the flocking dynamics\n",
"\n",
"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.\n",
"\n",
"### *Cohesion*: staying close to nearby flockmates\n",
"\n",
"

*Credit: Craig Reynolds (public domain)*\n",
"\n",
"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:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "b6488985",
"metadata": {},
"outputs": [],
"source": [
"def cohesion_force(positions, velocities, cohesion_strength=0.01):\n",
" return cohesion_strength * (positions.mean(axis=0)[np.newaxis] - positions)"
]
},
{
"cell_type": "markdown",
"id": "be117c86",
"metadata": {},
"source": [
"Now that we have defined our first force function, we can create a new animation with this force applied:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "6958bb69",
"metadata": {},
"outputs": [],
"source": [
"positions, velocities = initialise_boid_states(rng)\n",
"animate_flock(positions, velocities, [cohesion_force])"
]
},
{
"cell_type": "markdown",
"id": "b1de4b68",
"metadata": {},
"source": [
"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.\n",
"\n",
"### *Separation*: avoiding collisions with nearby flockmates\n",
"\n",
"

*Credit: Craig Reynolds (public domain)*\n",
"\n",
"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](https://numpy.org/doc/stable/reference/generated/numpy.where.html). Calling \n",
"```Python\n",
"values = np.where(conditions, values_if_true, values_if_false)\n",
"```\n",
"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)."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "e8391017",
"metadata": {},
"outputs": [],
"source": [
"def separation_force(positions, velocities, separation_strength=1., separation_distance=10.):\n",
" displacements = positions[np.newaxis] - positions[:, np.newaxis]\n",
" are_close = (displacements**2).sum(-1)**0.5 <= separation_distance\n",
" return separation_strength * np.where(are_close[..., None], displacements, 0).sum(0)"
]
},
{
"cell_type": "markdown",
"id": "e8441379",
"metadata": {},
"source": [
"We can now run and animate a simulation with both the forces we have defined:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "c4ba737f",
"metadata": {},
"outputs": [],
"source": [
"positions, velocities = initialise_boid_states(rng)\n",
"animate_flock(positions, velocities, [cohesion_force, separation_force])"
]
},
{
"cell_type": "markdown",
"id": "a33596d6",
"metadata": {},
"source": [
"### *Alignment*: matching velocity with nearby flockmates\n",
"\n",
"

*Credit: Craig Reynolds (public domain)*\n",
"\n",
"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`:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "f827c944",
"metadata": {},
"outputs": [],
"source": [
"def alignment_force(positions, velocities, alignment_strength=0.125, alignment_distance=100):\n",
" displacements = positions[np.newaxis] - positions[:, np.newaxis]\n",
" velocity_differences = velocities[np.newaxis] - velocities[:, np.newaxis]\n",
" are_close = (displacements**2).sum(-1)**0.5 <= alignment_distance\n",
" return -alignment_strength * np.where(are_close[..., None], velocity_differences, 0).mean(0)"
]
},
{
"cell_type": "markdown",
"id": "98b49c13",
"metadata": {},
"source": [
"Now we are finally ready to run and animate a simulation with forces for all three rules defined:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "72186e72",
"metadata": {},
"outputs": [],
"source": [
"positions, velocities = initialise_boid_states(rng)\n",
"animate_flock(positions, velocities, [cohesion_force, separation_force, alignment_force])"
]
},
{
"cell_type": "markdown",
"id": "6b1f7faa",
"metadata": {},
"source": [
"## Conclusion and extensions\n",
"\n",
"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.\n",
"\n",
"If you are interested in exploring this model further here are some ideas for extensions\n",
"\n",
" * 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?\n",
" * 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?\n",
" * 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](https://matplotlib.org/stable/tutorials/toolkits/mplot3d.html) useful for the visualisation side.\n",
"\n"
]
}
],
"metadata": {
"jekyll": {
"display_name": "The Boids"
},
"jupytext": {
"main_language": "python",
"notebook_metadata_filter": "-kernelspec,jupytext,jekyll"
}
},
"nbformat": 4,
"nbformat_minor": 5
}