Just before I started my paternity leave I saw a documentary about flocking animals and how these can be modeled in computer simulation. I found it very exciting and eagerly thought I could implement a simulation in a days work. I, however, quite underestimated the free time available while on paternity leave and the project extended in to many pieces spread over a few weeks. First a little teaser from the animation I made using python and Matplotlib, where the Boids (birds, fish or flying sheep ..) are in blue and predators are in red:

The definition of the term Boids and the rules that describes the motion of flock animals are excellently described on Wikipedia:

Boids is an artificial life program, developed by Craig Reynolds in 1986, which simulates the flocking behaviour of birds. His paper on this topic was published in 1987 in the proceedings of the ACM SIGGRAPH conference. The name refers to a “bird-like object”, but its pronunciation evokes that of “bird” in a stereotypical New York accent.

As with most artificial life simulations, Boids is an example of emergent behavior; that is, the complexity of Boids arises from the interaction of individual agents (the boids, in this case) adhering to a set of simple rules. The rules applied in the simplest Boids world are as follows:

Separation: steer to avoid crowding local flockmates

Alignment: steer towards the average heading of local flockmates

Cohesion: steer to move toward the average position (center of mass) of local flockmates

Below a youtube video of a natural phenomenon called “Sort Sol”, Black sun, where a very large group of common starlings blackens the skies over Ribe in Denmark:

## Animation using matplotlib and ffmpeg

The following video is the result of a run with 960 timesteps, 400 boids and 5 predators. I have used several extra rules and the predators in this run are dumb in the sense that they move in a predefined pattern (the Trefoil Knot for those interested), but it is quite easy to make them hunt the boids as the code is there, but the weight for this rule is 0. The behavior of the boids are controlled by the weights of the rules relative to each other.

## Python Implementaion

I have implemented the rules with animation in python using matplotlib and ffmpeg on Linux. The source is printed here, but can also be downloaded directly:

boids-animation.py 8.0KB.

#!/usr/bin/python import numpy import time import sys import math import matplotlib.pyplot as plt import matplotlib.animation as animation from mpl_toolkits.mplot3d import Axes3D ################# # Runtime parameters NumberOfBoids = 400 NumberOfPredators = 5 Dimensions = 3 MinSeparation = 7 framespersecond = 24 timesteps = 10 #Should be divisible by framespersecond # Predator settings PredatorRadius = 25 PredatorSight = 10 WeightCenterOfMassP = 0.00 WeightAttackBoid = 0 #2 WeightRandomP = 0.5 #1 WeightKnot = 2 # Boids behavior - weight on rules WeightCenterOfMass = 0.03 WeightSeparation = 1 WeightAlignment = 0.125 WeightRandom = 0.0 WeightAvoidPredator = 1 WeightCenter = 0.02 #0.03 MaxVelocityP = 1 MaxVelocity = 1 center = 50*numpy.ones(Dimensions) ########### class boid: """This class defines the boid and the rules of behavoir""" def __init__(self, dimension, index): self.index = index self.dim = dimension self.pos = numpy.random.uniform(0,100,self.dim) self.vel = numpy.zeros(self.dim) self.size = numpy.random.rand(1) def limitvelocity(self,maxvel): """Limiting the speed to avoid unphysical jerks in motion""" if numpy.linalg.norm(self.vel) > maxvel: self.vel = (self.vel/numpy.linalg.norm(self.vel))*maxvel def normalize(vector): """Normalizes a vector""" vector = vector/numpy.linalg.norm(vector) return vector def centerofpos(flock): """Calculates the 'center of mass' for the flock of boids""" com = numpy.zeros(Dimensions) for dim in range(Dimensions): for boid in flock: com[dim] = com[dim] + boid.pos[dim] com[dim] = com[dim]/len(flock) return com def centerofvel(flock): """Calculates the mean velocity vector for the flock of boids""" cov = numpy.zeros(Dimensions) for dim in range(Dimensions): for boid in flock: cov[dim] = cov[dim] + boid.vel[dim] cov[dim] = cov[dim]/len(flock) return cov def plot_raw(): """Sets up the basic plotting paramaters""" global ax ax = fig.add_subplot(111, projection='3d') ax.set_autoscale_on(False) ax.set_xlabel('X Axis') ax.set_ylabel('Y Axis') ax.set_zlabel('Z Axis') ax.set_xlim3d([0.0, 100.0]) ax.set_ylim3d([0.0, 100.0]) ax.set_zlim3d([0.0, 100.0]) def plot_init(): """The initial draw""" plot_raw() com = centerofpos(flock) boidpos = numpy.zeros((NumberOfBoids, Dimensions)) predpos = numpy.zeros((NumberOfPredators, Dimensions)) for index,boid in enumerate(flock): boidpos[index]= boid.pos bx = boidpos[:,0] by = boidpos[:,1] bz = boidpos[:,2] for index,predator in enumerate(predators): predpos[index]= predator.pos px = predpos[:,0] py = predpos[:,1] pz = predpos[:,2] def plot_update(timestep): """ This plot updates the plot with the new positions of the boids and predators""" fig.clf() plot_raw() # Execute a calculation of the next set of moves mainloop(timestep) # Plots the boids in space using matplotlib # http://matplotlib.org/examples/mplot3d/scatter3d_demo.html boidpos = numpy.zeros((NumberOfBoids, Dimensions)) predpos = numpy.zeros((NumberOfPredators, Dimensions)) # Plot center of mass as a black dot com = centerofpos(flock) ax.scatter(com[0], com[1], com[2], color='black') # Plot boids as blue dots for index,boid in enumerate(flock): boidpos[index]= boid.pos bx = boidpos[:,0] by = boidpos[:,1] bz = boidpos[:,2] ax.scatter(bx, by, bz, c='b') # Predators are red dots if len(predators) > 0: for index,predator in enumerate(predators): predpos[index]= predator.pos px = predpos[:,0] py = predpos[:,1] pz = predpos[:,2] ax.scatter(px, py, pz, s=100,color='red') def mainloop(timestep): """This is the calculation timeloop""" print "Timestep "+str(timestep) if len(predators) > 0: for predator in predators: # Predator RULE 1. Cohesion - Steer to move towoards the center of mass prule1 = numpy.zeros(Dimensions) prule1 = (centerofpos(flock) - predator.pos)*WeightCenterOfMassP # Predator RULE 2. (unused and now just a placeholder prule2 = numpy.zeros(Dimensions) # Predator RULE 3. Attack boids within range prule3 = numpy.zeros(Dimensions) for boid in flock: difference = predator.pos - boid.pos distance = numpy.linalg.norm(difference) if distance < PredatorSight: prule3 = prule3 - normalize(difference) prule3 = prule3*WeightAttackBoid # Predator RULE 4. Move the predator around in smooth way around the center of the cube # http://en.wikipedia.org/wiki/Trefoil_knot prule4 = numpy.zeros(Dimensions) t = (timestep/float(timesteps))*4*math.pi + predator.index*(math.pi/4.0) if Dimensions == 3: prule4[0] = (2.0+math.cos(3.0*t))*math.cos(2.0*t) prule4[1] = (2.0+math.cos(3.0*t))*math.sin(2.0*t) prule4[2] = math.sin(3*t) else: prule4 = numpy.zeros(Dimensions) prule4 = prule4 * WeightKnot # Move the predator predator.vel = prule1 + prule2 + prule3 + prule4 predator.limitvelocity(MaxVelocityP) predator.pos = predator.pos + predator.vel for boid in flock: # Boid RULE 1. Cohesion - Steer to move towoards the center of mass rule1 = numpy.zeros(Dimensions) rule1 = (centerofpos(flock) - boid.pos)*WeightCenterOfMass # Boid RULE RULE 2. Separation - steer to avoid crowding local flockmates rule2 = numpy.zeros(Dimensions) for boid2 in flock: difference = boid2.pos - boid.pos distance = numpy.linalg.norm(difference) if distance < MinSeparation and boid2 != boid: rule2 = rule2 - normalize(difference)/distance rule2 = rule2*WeightSeparation # Boid RULE 3. Alignment - Steer towards the average heading of local flockmates rule3 = numpy.zeros(Dimensions) rule3 = (centerofvel(flock) - boid.vel)*WeightAlignment # The following rules are custom rules just added for fun. # Boid RULE 4. Try to move towards the center of the grid rule4 = (center - boid.pos)*WeightCenter # Boid RULE 5. Add some randomness rule5 = numpy.random.uniform(-1,1,Dimensions)*WeightRandom # Boid RULE 6. Avoid the predator if len(predators) > 0: rule6 = numpy.zeros(Dimensions) for predator in predators: difference = predator.pos - boid.pos distance = numpy.linalg.norm(difference) if distance < PredatorRadius: rule6 = (rule6 - difference)*WeightAvoidPredator else: rule6 = numpy.zeros(Dimensions) # Move the boids boid.vel = rule1 + rule2 + rule3 + rule4 + rule5 + rule6 boid.limitvelocity(MaxVelocity) boid.pos = boid.pos + boid.vel #Generate the the flock of boids flock = [boid(Dimensions,count) for count in range(NumberOfBoids)] if NumberOfPredators == 0: predators = [] else: predators = [boid(Dimensions,count) for count in range(NumberOfPredators)] # HACK - centered starting positions for predator in predators: predator.pos = numpy.random.uniform(45,55,Dimensions) # Define the figure: fig = plt.figure() plot_raw() anim = animation.FuncAnimation(fig, plot_update, init_func=plot_init,frames=timesteps, interval=20, blit=True) anim.save('basic_animation.mp4', fps=framespersecond, codec='mpeg4', clear_temp=True, frame_prefix='_tmp') |

## References and further reading

- Conrad Parker, a PhD student in Computer Science from the University of Kyoto has a great guide on how to implement the rules in pseudocode:

http://www.vergenet.net/~conrad/boids/pseudocode.html - The Blender project has interesting points on further rules that can be implemented:

http://www.blender.org/development/release-logs/blender-246/particles/boids-physics/ - The documentation for the animation part of the matplotlib

http://matplotlib.org/examples/animation/index.html

http://matplotlib.org/api/animation_api.html

Did you choose the rules’ weights empirically??

The weights were just estimate to make it look reasonable, but I imagine they depend a bit on the setup, so do not think they very universal.