F# Flocking (boids)

A couple of years ago I was blogging about some simple boid behaviour I had implemented in C# (here). I came across this paper which had some interesting ideas I wanted to try out: "Self-organised complex aerial displays of thousands of starlings: a model" by H. Hildenbrandt, C. Carere, C-K. Hemelrijk.

As I’m getting more into F#, I thought this would be a simple reasonable sized application to create to get a feel for how F# works on a slightly larger project (though definitely not large-scale).

Before taking a look at the code, I’ll quickly describe which parts of the paper I implemented.

The video generated is available on http://www.youtube.com/watch?v=eil5K7Ir3i8

I haven’t implemented all of the ideas in the paper, so my simple application doesn’t exhibit the same realistic flocking behaviour, but it does have more realistic behaviour than my earlier efforts. The interesting behaviours were the desire to not stray too far from the roost, the attempt to maintain a steady cruise speed, and to use the nearest seven topological neighbours for cohesion and alignment. I implemented the simulation in 2D for simplicity.

I hardcoded a perception radius instead of implementing the continual adjustment of the perception radius for seven topological neighbours. I also have totally omitted the flight dynamic simulation (no gravity, no banked turns). The paper discusses that the banked turns gives a further air of realism to the simulation.  

The coding

This was a real pleasure to implement in F# – a lot less time was thinking about application design, classes and their interactions, the logic in the paper was easy to transfer directly to the keyboard. F# is very compact, so the core logic can be seen on just a couple of screens.

I didn’t feel the lack of intellisense too much when developing the core algorithm, I did miss the ability to refactor functions and navigate to usages. I especially missed intellisense when developing the WPF visualisation part of the app, when interacting with .NET objects; I did miss the ability for Visual Studio to automatically add the required namespaces. I must have being spoilt by Visual Studio and Resharper for too long!

The actual WPF application wasn’t such a good fit for F# – there’s no generated code behind file for the XAML, and I feel that using F# would be painful in a WPF or Silverlight application (but just for the view, it should be OK for the ViewModels down).

I implemented a F# Vector type, which can be specified with units of measure (or none). I used units of measure throughout – this was really powerful, and did let me quickly find a few bugs in the implementation.

(NOTE: I still haven’t found a code websnippet tool I’m happy with – you need to click in each of the regions below and scroll down and right to see the whole code). Alternately, the zipped up solution can be downloaded from http://www.taumuon.co.uk/jabuka/FSharpFlock.zip 


module Vector3D

type Vector3D<[<Measure>] 'u>(x : float<'u>, y : float<'u>, z : float<'u>) =
static member Zero() = Vector3D<_>(0.0<_>, 0.0<_>, 0.0<_>)

member v.X = x
member v.Y = y
member v.Z = z

static member (+) (lhs:Vector3D<_>, rhs:Vector3D<_>) =
Vector3D(lhs.X + rhs.X, lhs.Y + rhs.Y, lhs.Z + rhs.Z)

static member (-) (lhs:Vector3D<_>, rhs:Vector3D<_>) =
Vector3D(lhs.X - rhs.X, lhs.Y - rhs.Y, lhs.Z - rhs.Z)

static member (*) (v:Vector3D<_>, a:float<_>) =
Vector3D(v.X * a, v.Y * a, v.Z * a)

static member (*) (a:float<_>, v:Vector3D<_>) =
Vector3D(a * v.X, a * v.Y, a * v.Z)

static member (/) (v:Vector3D<_>, a) =
Vector3D(v.X / a, v.Y / a, v.Z / a)

member v.DotProduct(rhs:Vector3D<_>) = (v.X * rhs.X) + (v.Y * rhs.Y) + (v.Z * rhs.Z)

member v.magnitude = sqrt(v.DotProduct(v)) * 1.0<_>

member lhs.CrossProduct(rhs:Vector3D<_>) =
Vector3D((lhs.Y * rhs.Z - lhs.Z * rhs.Y) * 1.0<_>,
(-lhs.X * rhs.Z + lhs.Z * rhs.X) * 1.0<_>,
(lhs.X * rhs.Y - lhs.Y * rhs.X) * 1.0<_>)

member v.normalise =
let magnitude = float v.magnitude
Vector3D<_>((v.X / magnitude), (v.Y / magnitude), (v.Z / magnitude))

let sumVectors(vectors : Vector3D<_>[]) =
let initial = Vector3D<_>(0.0<_>, 0.0<_>, 0.0<_>)
Array.fold (+) initial vectors


module BoidUtils

open Microsoft.FSharp.Math
open Vector3D
open SI

let radiusRoost = 150.0<m>
let hardRadius = 2.0<m> // 0.2<m>
let mass = 0.08<kg>
let timeStep = 0.005<s>
let relaxationTime = 0.05<s>
let cruiseSpeed = 20.0<m/s>
let horizontalRoostWeighting = 0.01<N/m>
let weightingAlignment = 0.5<kg * s^-2>
let weightingCohesion = 1.0<kg s^-2>
let weightingSeparation = 2.0
let perceptionRadius = 50.0<m>

type BodyAxes =
{ Forward:Vector3D<1>;
Up:Vector3D<1> }

type Boid =
{ Position:Vector3D<m>;
Speed:float<m / s>;
Orientation:BodyAxes; }

// All parameterless functions are evaluated once, just on module opening, so pass random in.
let InitialiseRandomPosition(rand:System.Random) =
Vector3D<m>((300.0 * (rand.NextDouble()-0.5)) * 1.0<m>,
((300.0 * (rand.NextDouble()-0.5)) * 1.0<m>),

let InitialiseRandomVelocity(rand:System.Random) =
Vector3D<m/s>((100.0 * (-0.5 + rand.NextDouble()) * 1.0<m/s>),
(100.0 * (-0.5 + rand.NextDouble())) * 1.0<m/s>,

let InitialiseRandomOrientation(rand:System.Random) =
{Forward=Vector3D(0.0, 1.0, 0.0);
Side=Vector3D(1.0, 0.0, 0.0);
Up=Vector3D(0.0, 0.0, 1.0)}

let setOrientation(oldOrientation:BodyAxes, velocity:Vector3D<m/s>) =
let normalisedVelocity = velocity.normalise
let y = normalisedVelocity.CrossProduct(Vector3D<m / s>(0.0<m/s>, 0.0<m/s>, 1.0<m/s>))
{oldOrientation with Forward=normalisedVelocity * 1.0<m^-1 s>; Side=y*1.0<m^-2 s^2>}

let calculateCruiseSpeedForce (boid:Boid) =
(mass / relaxationTime) * (cruiseSpeed - boid.Speed) * boid.Orientation.Forward

let calculateRoostForce (boid:Boid) =
let horizontalPosition = Vector3D(boid.Position.X, boid.Position.Y, 0.0<_>)
let distanceFromOrigin = horizontalPosition.magnitude
match (distanceFromOrigin) with
| _ when distanceFromOrigin < radiusRoost -> Vector3D<N>(0.0<N>, 0.0<N>, 0.0<N>)
| _ -> let normalRoostingArea = horizontalPosition.normalise
let d = boid.Orientation.Forward.DotProduct normalRoostingArea
let distanceFromRoost = distanceFromOrigin - radiusRoost
let orientationRoostDotProduct = boid.Orientation.Side.DotProduct normalRoostingArea
let weightingFactor = match (orientationRoostDotProduct) with
| _ when orientationRoostDotProduct > 0.0<m> -> -1.0
| _ -> 1.0
weightingFactor * (radiusRoost * horizontalRoostWeighting * (0.5 + (0.5<m^-1> * d)) * (boid.Orientation.Side))

let findDistanceBetweenBoids(boid:Boid, other:Boid) =
(boid.Position - other.Position).magnitude

let findNearestNeighbours(boid:Boid, boids:Boid[]) =
let sortedByDistance = boids |> Array.sortBy(fun other -> findDistanceBetweenBoids(boid, other))
Array.sub sortedByDistance 0 7

let findAverageForwardDirectionDifference(boid:Boid, boids:Boid[]) =
let differences = boids |> Array.map (fun i -> 1.0<m> * (i.Orientation.Forward - boid.Orientation.Forward))
let sumDifferences = sumVectors(differences)
(1.0 / (float sumDifferences.magnitude)) * sumDifferences

let calculateAlignmentForce(boid:Boid, nearest:Boid[]) =
let averageDifference = findAverageForwardDirectionDifference(boid, nearest)
weightingAlignment * averageDifference

let findAveragePosition(boid:Boid, boids:Boid[]) =
let positions = boids |> Array.map (fun i -> i.Position)
let sumPositions = sumVectors(positions)
(1.0 / float boids.Length) * sumPositions

let findNeighboursInRadius(boid:Boid, boids:Boid[], radius:float<m>) =
boids |> Array.filter(fun other -> other <> boid && findDistanceBetweenBoids(boid, other) <= radius)

let calculateCentrality(boid:Boid, boids:Boid[]) =
let separations = boids |> Array.map(fun i -> (i.Position - boid.Position).normalise)
let sumSeparations = sumVectors(separations)
let count = boids.Length
match (count) with
| 0 -> 1.0
| _ -> (1.0 / float count) * (sumSeparations.magnitude / 1.0<m>)

let calculateCohesionForce(boid:Boid, nearest:Boid[], boidsInPerceptionRadius:Boid[]) =
let boidsOutsideHardRadius = nearest |> Array.filter(fun i -> abs ((boid.Position - i.Position).magnitude) > hardRadius)
let centrality = calculateCentrality(boid, boidsInPerceptionRadius)
let averagePosition = findAveragePosition(boid, nearest)
centrality * weightingCohesion * (averagePosition - boid.Position)

let calculateSeparationForce(boid:Boid, boidsInPerceptionRadius:Boid[]) =
let nearest = boidsInPerceptionRadius
let separations = nearest |> Array.map(fun i -> i.Position - boid.Position)
let sigma = 1.8
let forcesToNeighbours = separations |> Array.map(fun i ->
let magnitude = i.magnitude
let multiplier =
match (magnitude) with
| _ when magnitude < hardRadius -> 1.0
| _ -> System.Math.Exp(-((magnitude - hardRadius)*(magnitude - hardRadius)/1.0<m^2>) / (sigma * sigma))
multiplier * magnitude * (i.normalise) * 1.0<kg * m^-1 * s^-2>)
let sumForces = sumVectors(forcesToNeighbours)
match (nearest.Length) with
| _ when (nearest.Length) = 0 -> Vector3D<N>.Zero()
| _ -> (-weightingSeparation / float nearest.Length) * sumForces

let calculateSocialForce(boid:Boid, boids:Boid[]) =
let nearest = findNearestNeighbours(boid, boids)
let boidsInPerceptionRadius = findNeighboursInRadius(boid, boids, perceptionRadius)
calculateAlignmentForce(boid, nearest)
+ calculateCohesionForce(boid, nearest, boidsInPerceptionRadius)
+ calculateSeparationForce(boid, boidsInPerceptionRadius)

let calculateForce (boid:Boid, boids:Boid[]) =
(boid |> calculateRoostForce)
+ (boid |> calculateCruiseSpeedForce)
+ (calculateSocialForce(boid, boids))

let iterateBoid (boid:Boid, boids:Boid[]) =
let originalPosition = boid.Position
let originalVelocity = boid.Speed * boid.Orientation.Forward
let force = calculateForce(boid, boids)
let acceleration = force/mass
let velocity = originalVelocity + (acceleration * timeStep)
let position = originalPosition + (velocity * timeStep)
let newOrientation = setOrientation(boid.Orientation, velocity)



module Boids

open SI
open Vector3D
open BoidUtils

open System
open System.IO
open System.Windows
open System.Windows.Threading
open System.Windows.Controls
open System.Windows.Shapes
open System.Windows.Media
open System.Windows.Media.Imaging

let window = Application.LoadComponent(new System.Uri("/FSharpFlock;component/MainWindow.xaml",
System.UriKind.Relative)) :?> Window

let rand = new System.Random()
let viewWidth = 480.0
let viewHeight = 360.0

let mutable frameCount = 0
let mutable boids =
[|for i in 0 .. 300 ->
let position = InitialiseRandomPosition(rand)
let velocity = InitialiseRandomVelocity(rand)
let tempOrientation = {Forward=Vector3D(0.0, 1.0, 0.0);
Side=Vector3D(1.0, 0.0, 0.0);
Up=Vector3D(0.0, 0.0, 1.0)}
let orientation = setOrientation(tempOrientation, velocity)
{Position = position; Speed=velocity.magnitude; Orientation=orientation}

let updateBoids(boids) =
boids |> Array.map (fun boid -> iterateBoid(boid, boids))

let (?) (fe:FrameworkElement) name : 'T =
fe.FindName(name) :?> 'T

let GetRotationForBoid(boid:Boid) =
let forward = boid.Orientation.Forward
let angleRadians = Math.Atan2(forward.X, forward.Y)
let angleDegrees = angleRadians * 180.0 / Math.PI
let rotateTransform = new Media.RotateTransform(angleDegrees, 0.0, 0.0)

let color = Media.Color.FromArgb(64uy, 255uy, 0uy, 0uy)

let saveFrame(canvas : Canvas) =
let size = new Size(viewWidth, viewHeight)
canvas.Arrange(new Rect(size))
let renderTargetBitmap = new RenderTargetBitmap(int viewWidth, int viewHeight, 96.0, 96.0, PixelFormats.Pbgra32)
let sourceBrush = new VisualBrush(canvas)

let bitmapFrame = BitmapFrame.Create(renderTargetBitmap)
let jpegEncoder = new JpegBitmapEncoder()
let filename = System.String.Format("C:Anim{0:0000}.jpg", frameCount)
use stream = new FileStream(filename, FileMode.CreateNew)

let createBoidGraphics(boid:Boid) =
let obj : Polygon = new Polygon()
obj.Points <- new Media.PointCollection( [|new Point(-10.0, -5.0);
new Point(0.0, 0.0);
new Point(-10.0, 5.0);
new Point(-5.0, 0.0)|] )
obj.RenderTransform <- GetRotationForBoid(boid)
obj.Fill <- new Media.SolidColorBrush(color)
obj.Stroke <- Media.Brushes.Black
obj.StrokeThickness <- 1.0

let drawBoids() =
let win : Window = window
let canvas : Canvas = win?canvas
for i = 0 to boids.Length - 1 do
let graphicalBoid = createBoidGraphics(boids.[i])
let unitlessPosition = boids.[i].Position * 1.0<m^-1>
System.Windows.Controls.Canvas.SetTop(graphicalBoid, (viewHeight / 2.0) + (unitlessPosition.X))
System.Windows.Controls.Canvas.SetLeft(graphicalBoid, (viewWidth / 2.0) + (unitlessPosition.Y))
canvas.Children.Add(graphicalBoid) |> ignore
// saveFrame(canvas)
frameCount <- frameCount + 1

let timer = new DispatcherTimer();
timer.Tick.Add(fun evArgs ->
|> Seq.iter(fun x -> boids <- boids |> updateBoids)
timer.Interval <- TimeSpan.FromMilliseconds(1.0);

let setup(win: Window) =
win.Loaded.Add(fun _ -> timer.Start())

let main args =
setup window
(new Application()).Run(window) |> ignore


<Window xmlns="http://schemas.microsoft.com/winfx/2006/xaml/presentation"
Title="F# Flock" Height="360" Width="480">
<Canvas Name="canvas" HorizontalAlignment="Stretch" VerticalAlignment="Stretch"

The rest of the code I hope is simple to read and understand. Coding this up I had some of the most fun programming for a while. With F#, I found I was able to concentrate on the core of the algorithm and very quickly see results.

For this type of project, very algorithmic, F# was a perfect fit. I’d like to get a feel for how to implement F# on a large scale project, and to see how it feels to expose the functionality via an object oriented layer.

Faster boids

I used the built in performance profiler in VSTS to look at speeding the boids up (and there seems to be a bug in the profiler – when memory profiling was turned on, it seemed unable to resolve symbols, so I switched to the CLRProfiler to look at the memory being used).

Most of the time in the application was spent in calculating the boids’ steering behaviours.
The optimisations here were to try to reduce the time spent in calculating the forces – one optimisation was that instead of each behaviour (cohesion, separation, alignment) each scanning through the entire set of boids, the results of one of the behaviours could be used by both others.

Secondly, the cohesion force was not updated on each frame – the value is cached for 10 iterations. This should not be a fixed value – there should be some heuristic to determine how it should be varied once work is in to measure the “quality” of the flock.

This boosted the frame rate from approx 3fps to 10fps (dropping to 5 fps after some time).

Next came some micro-optimisations.

Changing the Vector to be a struct instead of a class (to avoid over-stressing the garbage collector) resulted in 14fps dropping to 8 fps after running.

Other micro-optimisations that I tried, but didn’t result in any real performance difference included:

Making the trail generator to use a fixed size array (instead of adding/removing from a list).

Removing IEulerRigidBody, calling onto objects directly instead of through an interface.

Changing all doubles to be floats (didn’t actually keep this change in the code).

For Scene8, the drawing code was the real bottleneck, changing the gluSphere to be drawn using VertexBufferObjects resulted in double the frame rate (and it isn’t particularly optimised – could use glDrawElements).

For Scene7, it’s limited by the O(n2) collision detection (this is the next thing to tackle).


An improvement of the flocking behaviour in Jabuka is available.

There are now very simplistic birds drawn, to give an idea of the boid’s orientation.

The separation behaviour has been fixed. The behaviour was garbage before, the force away from another boid was proportional to the distance away from it. Now the force increases the nearer another boid is.

The boid’s orientation is taken as the smallest angle between being aligned with the World’s x axis, and the boid’s velocity. The roll doesn’t behave realistically (see scene 13), but the boid’s behaviour isn’t realistic. A boid would really change its velocity by adjusting its orientation (change the wing shape to apply a torque to allow it to rotate…)

This paper describes how the steering behaviours should change the intention of the boid, instead of providing actual forces onto the boid.

A simple particle system has been introduced to allow the path of a few individual boids to be more easily seen.

Running under the Visual Studio profiler, most of the time was seen to be spent finding the flockmates in proximity of the boid.

The same calculation was being done twice, for both the cohesion and alignment behaviours. Adding a ProximityFinderCache increased the frame rate by 50%, with no change in the calculations. Adjusting the code so that the boids in proximity are only calculated on every third iteration added another 20% to the framerate (giving 15fps with only two TrailGenerators attached).