tl;dr – code now availalbe on github

I previously posted an optimized version of Craig Reynolds flocking algorithm. Nothing very fancy – simply adding a way to store the boid objects in a grid data structure which makes looking up the distances between boids much faster. I felt the need to do this because there are plenty of boid examples out there in various languages/frameworks and at the end of most of them they mention how it could be greatly optimized etc etc but never actually give any examples on how to do it.

So I created an example in java using pulpcore. Recently however i’ve been playing with python, and I love it. I also realize that the java version is rather complicated and so doesn’t work too well as an example for others. Therefore i’ve re-implemented the algorithm in python, using pyglet for display. Pyglet is a wrapper for opengl written in pure python. Its pretty low-level but that’s not important – the optimized boid simulation is very much separated from the graphics part.

The simulation has 3 important classes. The Boid class is for the boids themselves, the BoidSwarm which is the grid data structure for holding the boids, and the FlockSimulation which just ties the two together and makes the data structure holding the boids completely separate from the boids themselves. In fact you could replace the BoidSwarm with a simple array and the simulation would still work, though more slowly.

Note also there are some important differences compared to the Java version because of the nature of Python. For example there’s no use of linked lists, since the builtin python list is faster than any alternatives. Similarly the buitin python hashmap is used rather than the lower level constructs of java.

[UPDATE] I did some very approximate benchmarking comparing the optimized data structure to a simple list approach. With 100boids on a 2010 13″ macbook I get around 45FPS optimized vs 10FPS for the unoptimized version. This is nice to see, since it makes the simulation usable in python, even though that’s not normally the primary choice for this kind of thing. I’m also sure that could be further improved. The main thing really is just to illustrate the algorithm for those who want it.

Anyway enough talking: some code. The files are available on dropbox here. This includes engine.py, which is the pyglet based graphics app, optboids.py which has the simulation code, and euclid.py (from cocos2d) which includes a required 2D vector class. Pyglet should be very easy to install since it has no dependencies. One caveat however: you might have issues getting it to work with virtualenv on MacOS, it’s probably easier not to use them.

To launch, just do python engine.py in the folder containing the files.

Below, the key classes of the simulation:

Boid – the boid knows its own dynamics (how it interacts with other boids), it just needs a list of its neighbours given to its interact function, which calculates the boids current acceleration vector. The update function takes a time delta (in seconds) and uses the acceleration to compute the velocity and new position.

Note that many Boid examples have separate functions for alignment, separation, and cohesion, however this adds a bit of unnecessary overhead which is avoided by melding the three into the interact function.

The trickiest part is actually getting all the parameters right – steering happens by adding vectors from the 3 effects, getting a good balance is a matter of fiddling with the factors.

BoidSwarm – the data structure. Includes a list which holds all the Boids, and a dict. The keys of this dict are tuples with the integer index of each grid square, and the values are lists which contain the boids in that square. There are also a few util functions for inserting a point in the correct grid square, as well as getting all the boids within a given radius around a point.

FlockSimulation – this just ties Boid and BoidSwarm together with an update function. It runs through the list of boids and passes the set of neighbouring boids to each boids interact function.

def limit(vector,lim): """ limit a vector to a given magnitude this is an 'in place' function, modifies the vector supplied. """ if abs(vector) > lim: vector.normalize() vector *= lim class Boid(object): """ Boids class """ influence_range = 40 minsep = 25.0 damp = 90.0 max_force = 20.0 max_speed = 100.0 max_steering_force = 20.0 drag = 0.9 cohesion_strength = 1.0 align_strength = 1.0 sep_strength = 1.0 cohesion_strength *= max_steering_force align_strength *= max_steering_force sep_strength *= max_steering_force # Get and set the position as a vector def _get_pos(self): return Vector2(self.x,self.y) def _set_pos(self,p): self.x = p.x self.y = p.y position = property(_get_pos,_set_pos) # Get and set the speed as a scalar def _get_speed(self): return abs(self.velocity) def _set_speed(self,s): if abs(self.velocity) == 0: velocity = Vector2(1,0) self.velocity.normalize() self.velocity *= s speed = property(_get_speed,_set_speed) #get and set the rotation as an angle def _get_rotation(self): return atan2(self.velocity.y,self.velocity.x) def _set_rotation(self,r): old_speed = self.speed #set the direction as a unit vector velocity.x = cos(r) velocity.y = sin(r) self.speed=old_speed rotation = property(_get_rotation,_set_rotation) def __init__(self,x, y): """ create a new boid at x,y """ self.x = x self.y = y self.acceleration = Vector2(0,0) self.velocity = Vector2(random.uniform(-self.max_speed, self.max_speed), random.uniform(self.max_speed, self.max_speed)) def __repr__(self): return 'id %d'%self.id def borders(self, top, bottom, left, right): """ Assume a square field """ if self.x < left: self.velocity += Vector2(1,0) * self.speed / 5 if self.y < top: self.velocity += Vector2(0,1) * self.speed / 5 if self.x > right: self.velocity -= Vector2(1,0) * self.speed / 5 if self.y > bottom: self.velocity -= Vector2(0,1) * self.speed / 5 def update(self,t): """ Method to update position by computing displacement from velocity and acceleration """ self.velocity += self.acceleration * t *self.drag if self.speed > self.max_speed: self.velocity *= 0.9 #limit(self.velocity,self.max_speed) self.position += self.velocity * t def interact(self, actors): """ Unit-unit interaction method, combining a separation force, and velocity alignment force, and a cohesion force """ separate_force = Vector2(0,0) align_force = Vector2(0,0) cohesion_sum = Vector2(0,0) cohesion_force = Vector2(0,0) count = 0 for other in actors: #vector pointing from neighbors to self diff = self.position - other.position d = abs(diff) #Only perform on "neighbor" actors, i.e. ones closer than arbitrary #dist or if the distance is not 0 (you are yourself) if d > 0 and d < self.influence_range: count = count+1 diff.normalize() if d < self.minsep: diff *= self.max_steering_force else: diff /= d # Weight by distance cohesion_sum += other.position # Add position separate_force += diff #Align - add the velocity of the neighbouring actors, then average align_force += other.velocity if count > 0: #calc the average of the separation vector separate_force /=count separate_force *= self.sep_strength #calc the average direction (normed avg velocity) align_force /=count align_force *= self.align_strength #calc the average position and calc steering vector towards it cohesion_sum /= count cohesion_force = self.steer(cohesion_sum,True) cohesion_force *= self.cohesion_strength #finally add the velocities sum = separate_force + cohesion_force + align_force #smooth the acceleration var self.acceleration = 0.2*self.acceleration + sum def steer(self, desired, slowdown=False): """ A helper method that calculates a steering vector towards a target If slowdown is true the steering force is reduced as it approaches the target """ desired -= self.position d = abs(desired) #If the distance is greater than 0, calc steering (otherwise return zero vector) if d > 0: desired.normalize() if slowdown and (d < 30.0): desired *= self.max_steering_force*d/self.damp else: desired *= self.max_steering_force steer = desired - self.velocity else: steer = Vector2(0,0) return steer def seek(self, target, m=1): self.acceleration += self.steer(target,False) self.acceleration = self.acceleration * m

class BoidSwarm(object): """ The grid to hold the boids """ def __init__(self,width, cell_w): """ Create data structure to hold the things. At the base is a table of cell nodes each containing an arbitrary number of units. Need whole number of cells, so cell width is the total width (for now in pixel units) divded by the divisions. The structure is always square. It can provide boundaries though good level design should mean you don't ever hit them. """ self.boids = [] #list of all the boids divs = int(floor(width/cell_w)) self.divisions = divs self.cell_width = cell_w self.cell_table = {} self.num_cells = divs*divs for i in range(divs): for j in range(divs): self.cell_table[(i,j)] = [] def add_boid(self,b): self.boids.append(b) self._insert(b) def _insert(self, b): c = self.find_cell_containing(b.x, b.y) c.append(b) def cell_num(self,x, y): """Forces units into border cells if they hit an edge""" i = int(floor(x / self.cell_width)) j = int(floor(y / self.cell_width)) #deal with boundary conditions if i < 0: i=0 if j = self.divisions: i = self.divisions-1 #zero indexing if j >= self.divisions: j = self.divisions-1 #zero indexing return (i,j) def find_cell_containing(self, x, y): """returns the cell containing a position x,y""" return self.cell_table[self.cell_num(x,y)] def find_near(self, x, y, influence_range): """return objects within radius influence_range of point x,y""" if influence_range == 0: _nearObjects = self.find_cell_containing(x,y) elif influence_range <= self.cell_width: _nearObjects = self.find_neighbour_cells(x,y) else: ext = ceil(influence_range/self.cell_width) _nearObjects = self.find_extended(x,y,ext) return _nearObjects def find_neighbour_cells(self, x, y): return self.cell_table[self.cell_num(x,y)] def find_extended(self, x, y, d): """ use to find set of cells surrounding the cell containing position x,y """ I,J = self.cell_num(x,y) d=int(d) group = [] for i in range(I-d,I+d): for j in range(J-d,J+d): if (i,j) in self.cell_table: group += self.cell_table[(i,j)] #list concat return group def rebuild(self): for cell in self.cell_table.values(): cell[:] = [] for b in self.boids: self._insert(b)

class FlockSimulation(object): """ Ties the BoidSwarm with the boids. boidswarm just holds the data, boids know how to interact with each other. This class keeps the two separated """ def __init__(self,starting_units = 100, field_size = 800): """ """ self.swarm = BoidSwarm(field_size+2*40,Boid.influence_range+5) self.field_size = field_size self.pad = self.swarm.cell_width #use to keep boids inside the play field for i in range(starting_units): b = Boid(random.uniform(100, 600), random.uniform(100, 600)) self.swarm.add_boid(b) self._cumltime = 0 #calculation var def update(self,dt): """dt is in seconds""" for b in self.swarm.boids: _close_boids = self.swarm.find_near(b.x,b.y,b.influence_range) b.interact(_close_boids) b.update(dt) w=self.field_size p=self.pad b.borders(p,w-p,p,w-p) #keep the boids inside the borders #rebuild the swarm once we've updated all the positions self.swarm.rebuild()

*Computers*

Nicolas Rougier

March 23, 2012

Nice article.

I’ve also coded ecently some boids using python and numpy (that give easy access to kdtrees). I will try to check if your solution can speedup the code (currently, I have around 1500 boids).

Code is available from http://code.google.com/p/glumpy/ (scroll down for screenshot/movie).

Nicolas