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

Vector3D.fs:

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

BoidUtils.fs:

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>;  
    Side: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>),   
             0.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>,  
                0.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)  
  {Position=position;Speed=velocity.magnitude;Orientation=newOrientation}

Program.fs:

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)  
  rotateTransform  
  
let color = Media.Color.FromArgb(64uy, 255uy, 0uy, 0uy)  
  
let saveFrame(canvas : Canvas) =   
  let size = new Size(viewWidth, viewHeight)  
  canvas.Measure(size)  
  canvas.Arrange(new Rect(size))  
  let renderTargetBitmap = new RenderTargetBitmap(int viewWidth, int viewHeight, 96.0, 96.0, PixelFormats.Pbgra32)  
  let sourceBrush = new VisualBrush(canvas)  
  renderTargetBitmap.Render(canvas)  
  
  let bitmapFrame = BitmapFrame.Create(renderTargetBitmap)  
  let jpegEncoder = new JpegBitmapEncoder()  
  jpegEncoder.Frames.Add(bitmapFrame)  
  let filename = System.String.Format("C:Anim{0:0000}.jpg", frameCount)  
  use stream = new FileStream(filename, FileMode.CreateNew)  
  jpegEncoder.Save(stream)  
    
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  
  obj  
  
let drawBoids() =  
  let win : Window = window  
  let canvas : Canvas = win?canvas  
  canvas.Children.Clear()  
  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 ->   
  [0..5]  
  |> Seq.iter(fun x -> boids <- boids |> updateBoids)  
  drawBoids()  
  )  
timer.Interval <- TimeSpan.FromMilliseconds(1.0);  
  
let setup(win: Window) =  
  win.Loaded.Add(fun _ -> timer.Start())  
  
[<STAThread>]  
[<EntryPoint>]  
let main args =  
  setup window  
  (new Application()).Run(window) |> ignore  
  0

MainWindow.xaml:

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

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.