Monte Carlo pricing of Exotic Options in F#

I’ve recently completed all exercises for the MOOC Monte Carlo methods in Finance, ran by the University of Madrid on Iversity. The course was excellent, both in its content and delivery.

The course started off by refreshing probability theory (cdf, pdf etc) and quickly moved onto generating random numbers, Brownian motion and stochastic differential equations followed by the pricing of options and calculating risk.
The final week’s homework is on the Monte Carlo calculation of a portfolio’s risk, using copulas with Student’s t distribution, but as that week’s homework is still running I’m blogging about the penultimate week’s homework.

This homework is about pricing path dependent options, and demonstrating how using a control variate can provide a better result for the same number of iterations (better meaning less error, i.e. with a smaller standard deviation). The homework further shows that when pricing a barrier
option, using an up-and-out option as the control variate produces better results than using a European option as the control variate.

The course delivery was extremely slick. The videos were filmed specially, this wasn’t a rehashing of a lecture in front of students (as seems to happen in some other online courses). The lectures were short, and to the point, managing to cram a lot of information in with zero fluff. I really appreciated this as I was watching on a daily commute.

The assessments were by multiple choice question, but with only one attempt, so it did end up actually feeling like a challenge. The questions were a mixture of maths and coding, with the balance moving more towards coding as the weeks went on.
The code was delivered in Matlab, which I had to get reacquainted to after a decade break. The main downside I found is that I ended up doing a fair bit of copy-paste coding, so completed some of the homeworks in .NET instead. I personally found the coding easy, and found the maths more challenging (in a good way).

Now, for the code – the F# code is available here: https://gist.github.com/taumuon/11302896

Apart from the Monte Carlo code, the fact that in F# all functions are in the curried form means that composing the payoff functions using partial function application is beautiful. This would be lambda hell to replicate in C#.

If this course runs again I’d wholeheartedly recommend it.

MOOCs

I’ve just enrolled on Iversity’s Monte Carlo Methods in Finance course, and have converted some of week 1’s Matlab demo code over to F# and Deedle: https://gist.github.com/taumuon/8602365

I spent the latter half of last year diving into various online courses. I completed Coursera’s Mathematical Methods for Quantitative Finance, and also followed along with a number of other courses to varying levels of completeness.

I completed three weeks assignments of Udacity’s Cuda programming course.  Week three was painful due to an error in the reference code, and week 4 was crashing due to a memory exception. I was using a GPU emulator in Ubuntu, and decided that it would be easier with real hardware. I watched the remaining videos and found the parallel algorithms explanations useful.

I completed the first two weeks of Coursera’s Scientific Computing. These were maths exercises, which I enjoyed and that’s what inspired me to do the Maths Methods for Quant Finance course. The Matlab exercises I was planning to do in F#, but left the course when other attendees were complaining that to complete the homework to the correct numerical accuracy you needed the exact version of Matlab the instructor was using, and they were unable to use Gnu Octave.

It is great that there’s so much free high standard material available. The fixed timescale nature of the course is a bit annoying – if work or life gets in the way one week it may make it impossible to catch up with the remainder of the course. I may get around to trying a course again next time it comes around though.

RunKeeper Visualisations

My last post described how I’m using the Twitter API to receive tweets off the live stream.

Since then, I’ve used the API to filter for tweets containing the #runkeeper hashtag, and used that to scrape the user’s activity from the RunKeeper site (including the GPS points from the user’s exercise). I’ve stored that information in a MongoDB, which has allowed me to do some simple visualisation:


The above video (best played at 720p) shows activities plotted against time of day (the sun overhead in the video indicates midday for that region).

I haven’t charted this to confirm, but to my eyes it looks like amount of exercise activity peaks around sunrise and sunset, with almost none at night time (which isn’t really a surprise).

For the curious, this is what the colour of the dots indicate:

let createSphere (latitude:float) (longitude:float) (activityType:string=
    let color = match activityType with
                | "RUN" -> "Red"
                | "WALK" -> "Green"
                | "BIKE" -> "Blue"
                | "MOUNTAINBIKE" -> "Orange"
                | "SKATE" -> "Yellow"
                | "ROWING" -> "Cyan"
                | "HIKE" -> "Brown"
                | "OTHER" -> "Black"
                | _ -> "Black"
    String.Format(spherelatitudelongitudecolor)

On a local scale, the actual GPS traces are of interest:

uk2

Activities are present in the more-populated regions of the UK. The blue and orange traces indicating cycling and mountain biking activities are clearly visible.

NOTE: if you think the map looks weird, it doesn’t have the usual aspect ratio – I’m not using a Mercator projection do plot the points, but am simply plotting the longitude and latitude linearly.

On an even smaller scale, landmarks around London are visible:

uk3

longoogle

The GPS data contains altitude information, so there are more interesting visualisations that could be done, including generating contour plots. Also, the above only contains three days worth of data – with a larger data set it would be possible to plot to determine whether activities peak on a weekend etc.

Implementation

Acknowledgements:

The 3D plot of the globe was drawn using POV-Ray. The 3D globe model was from grabcad, with the converted to a POV-Ray model using PoseRay.

The UK outline was obtained from the NOAA.

The code was implemented in F# (which was a pleasure to get back to after the C++ I’ve been doing recently). I did try the MongoDB.FSharp library to store my data records, but they failed to deserialise from the database. In any case, I wanted more control over the data types saved (I wanted to store the GPS data as GeoJson3DGeographicCoordinates, with the first point stored separately as a GeoJson2DGeographicCoordinates with an index on this value). I could have created .NET classes and used the BSON serializer, but it seemed more effort than writing directly to the DB (and this is about the only time I’ve seen the benefit of C# implicit conversions, but I can live without them in F#).

Why use the Twitter API, and why not scrape the RunKeeper site directly? That’s because of times – the RunKeeper website displays the activity time, but it is displayed in local time, and it’s not clear whether that’s in the user’s timezone, or the timezone of the activity. It seems cleaner to instead assume that the tweet has been posted as soon as the activity is finished and store that time as UTC (this assumption may of course not be true, but the results seem realistic).

Show me the Code!

The code’s not in a bad shape, but I would like to tidy it up a little before releasing it into the world. I’m busy with other things at the moment, but if I get much interest I can go ahead and do that…

Sipping from the twitter stream

I’ve been playing around with connecting to twitter’s streaming API, and displaying a live stream of tweets returned.
To do this, I was originally using DotNetOpenAuth, along with the HttpClient, which worked fine for the sample stream, but would return authentication errors for the filtered stream. I looked at the HTML message in fiddler2, and the oauth parameters weren’t ordered, which twitter requires. Instead, I’m using TwitterDoodle, which uses HttpClient.
The C#5/.NET4.5 async/await code is quite elegant – it’s not CPU intensive so it’s fine running on the UI thread, without blocking. My first instinct prior to C#5 would have been to use RX for this, but now if I’m doing something simple I’d stick to async/await, only using RX if doing something more complex like buffering or batching.

    private async void ProcessTweets()
    {
      using (var t = new TwitterConnector())
      {
        var response = await t.GetSampleFirehoseConnection();
        var res = await response.Content.ReadAsStreamAsync();
        using (var streamReader = new StreamReader(resEncoding.UTF8))
        {
          // streamReader.EndOfStream can block, so instead check for null
          while (!cts.IsCancellationRequested)
          {
            var r = await streamReader.ReadLineAsync();
            if (r == null) { return; }
            ProcessTweetText(r);
          }
        }
      }
    }

    private void ProcessTweetText(string r)
    {
      if (!string.IsNullOrEmpty(r))
      {
        var tweetJToken = JsonConvert.DeserializeObject<dynamic>(r);
        var tweetObj = tweetJToken["text"];
        if (tweetObj != null)
        {
          var tweetText = tweetObj.ToString();
          viewModel.Items.Add(tweetText);
        }
      }
    }

The equivalent F# async code obviously looks quite similar, with the added goodness of Type Providers. Time dependent, I am planning to do some more analysis of tweets which would be a good fit for F#.

type tweetProvider = JsonProvider<"SampleTweet.json"SampleList=true>

type MainWindowViewModel() =
  inherit ViewModelBase()
  let items = new ObservableCollection<string>()
  member x.Items
    with get () = items
  member x.ProcessTweet tweet =
    let tweetParsed = tweetProvider.Parse(tweet)
    match tweetParsed.Text with
    | Some(v->  x.Items.Add v
    | None -> ()
  member x.ProcessTweets =
    let a = async {
      use t = new TwitterConnector()
      let! response = t.GetSampleFirehoseConnection() |> Async.AwaitTask
      let! result = response.Content.ReadAsStreamAsync() |> Async.AwaitTask
      use streamReader = new StreamReader(resultEncoding.UTF8)
      let rec processTweetsAsync (s:StreamReader=
        async {
          let! r = s.ReadLineAsync() |> Async.AwaitTask
          // streamReader.EndOfStream can block, so instead check for null
          if r <> null then
            x.ProcessTweet r
            return! processTweetsAsync s
        }
      do! processTweetsAsync streamReader
    }
    a |> Async.StartImmediate
  member x.GoCommand = 
    new RelayCommand ((fun canExecute -> true), 
      (fun action -> x.ProcessTweets))

Option pricing in F# using the Explicit Finite Difference method

It’s been a little while since I’ve coded any F#, so I’ve done a little coding kata to polish off the rust.

I’ve converted the explicit finite difference option pricing example from (Paul Wilmott Introduces Quantitative Finance).

Call

I followed similar steps as those I performed in my previous blog post on the binomial tree option pricing; I converted the VBA code into very similar looking imperative F# code, before refactoring out the loops into recursive calls.

It was a useful exercise as converting the algorithm made me actually focus on the mechanics of how it worked, much more than simply translating VBA code into imperative F#. It’d be interesting to hear in the comments whether people find it easy to read.

The code is available https://gist.github.com/taumuon/4999749.

Binomial Tree Option Pricing in F#

This post will look at implementing the binomial tree algorithm in an imperative style, and refactoring it to be more functional.

The wikipedia page on the binomial options pricing model provides pseudo-code for the pricing of an American Put option.

That code can be implemented almost identically in F#.

clip_image002_thumb25255B125255D

and is called as follows:

clip_image004_thumb25255B125255D

Functional programmers have probably died a little on the inside on seeing that code, but one of the good thing about F# being multi-paradigm is that it’s easy to take code of this form and implement it, and test it for correctness before refactoring.

I did actually look at whether the results of the code were as expected at this point, but I’ll cover this later in this blog post.

The first thing to do is to change the initial values V to be created more functionally:

clip_image00625255B425255D

 

 

Now, the main algorithm consists of two parts, an outer loop where each ‘level’ of the tree is calculated stepping up towards the root, using the results of the previous level’s nodes (the nodes nearest to the leaves). The inner loop iterates over all of the nodes at a given level calculating the option’s value at that node.

First off, the outer loop can be replaced with a recursive method:

clip_image00825255B425255D

The next change is quite minor – it’s for each pass of the inner loop to create a new array, instead of mutating the vi array:

clip_image01025255B425255D

This change now means that there’s no need for any mutable state, so the initial values can be created as a sequence instead of an array, and the inner loop can be written in a functional style:

clip_image01225255B425255D

I find this code style is clearer to understand than following the array mutations to see what each part of the code is doing and how those parts interact.

If I was implementing this in C# I’d probably change the variable names to be more descriptive (e.g. strikePrice, riskFreeRate etc) but as one of the key advantages of F# is its succinctness it makes more sense to leave the well-known symbols in. Anyone implementing the algorithm with knowledge of the domain would know what the symbols mean without the more verbose names.

Validating the results

The calculated option value with the requested parameters is 4.486090. Running the EqutyOptions example in the .NET port of QuantLib, QLNet (http://sourceforge.net/projects/qlnet), with the same parameters gives the following results:

clip_image01425255B425255D

The figures are in the same ballpark, but the difference is curious.

Chapter 45 of GPU GEMS 2  creates the initial values in this way:

clip_image016_thumb25255B125255D

 

This gives a result of 4.486375 which still doesn’t match QLNet, but is nearer.

The result of running the Cox version of Bubo’s code from here http://cs.hubfs.net/topic/None/59053 also equals 4.486375, which matches the results of using this initialisation function. As an aside, I like how that examples shows how F# partial function application allows the variations on the algorithm to be composed, and much more succinctly than using OO.

Bubo’s Tian model does match QuantLib’s result exactly. I’ll have to dig a bit deeper to figure out which variations on the binomial tree model these different examples are using.

The code is available here: https://github.com/taumuon/FSharp-Binomial-Tree

More playing with monads

This blog post will look at using the continuation monad to sum the nodes in a tree, and then look at a case where a more complex monad is needed.

The book Real World Functional Programming gave an example of an unbalanced tree, which was deeply nested such that recursing over the tree to sum the values would blow the stack:

clip_image001_thumb25255B325255D

clip_image00225255B425255D

The solution provided was to sum the tree using continuation passing:

clip_image00325255B425255D

This is now able to sum the tree without crashing, but the code doesn’t flow as well as in the original sumTree method. However, monads can come to the rescue again, in this case the continuation monad:

clip_image00425255B425255D

The code flow now is easier to understand. The continuation passing style and continuation monad turns the code flow inside out, which isn’t a big deal, but means that the func to perform on the result needs to be passed in:

clip_image00525255B425255D

Parallel sum

This is all well and good, but to exploit multicore machines, I want to take advantage of the parallelism in the tree summation algorithm, by summing the different branches of the nodes using tasks. A naïve implementation to do this may be:

clip_image00625255B425255D

 

The problem with this is that it runs slower than the sequential case for a large balanced tree –if there are only 4 cores in the system it’s counter-productive to create thousands of tasks. The current depth of the tree should be passed around, and tasks spawned until traversal reaches that depth:

clip_image007_thumb25255B225255D

 

 

This is now faster for a large balanced tree, but for the deeply nested unbalanced tree it throws a StackOverflowException. Note that sumTreeParallel coped with arbitrarily nested trees, as the work to be done no longer kept being pushed onto the stack.

More state is being passed around (the current depth), which is external to the main business logic – it does feel like the state monad could be used here. First though, this is fixed up so that it has the efficiency of the parallel summation, but can cope with arbitrarily deeply nested trees:

clip_image00825255B425255D

Urghhhh. The code is definitely more convoluted now – I played with using the state monad to pass the depth around, but it’s not tail recursive. I thought that it might be possible to make use of the continuation monad here, but the problem is, is that it really needs aspects of both the state and continuation monads. I’m not sure how monad combinations work in F#, but I thought I’d throw this out here to see if anyone shows me the light. The post The Mother of All Monads does say that all other monads could be implemented in terms of the continuation monad, but I couldn’t see an example of anyone implementing the state monad in terms of the continuation monad.

State Monad in C# and F#

My colleagues and I recently have been trying to gain a deeper understanding of monads, other than using LINQ and RX. Watching all the channel9 videos on monads, the concept on the surface seemed quite simple, but we were felt that we were somehow missing some key point. After reading various blog posts, the penny is finally dropping.

Some post I read somewhere (I can’t remember where) said that it’s easier to understand the individual monads before trying to understand the whole generalised concept. Mike Hadlow’s blog series (http://mikehadlow.blogspot.co.uk/2011/01/monads-in-c1-introduction.html) really explain well the identity and maybe monads (which we already understood, which is just as well as we’re using the maybe monad in our codebase). Unfortunately he didn’t explain the state monad.

The only information we could find on the state monad in C# was Brian Beckman’s Channel 9 video: http://channel9.msdn.com/Shows/Going+Deep/Brian-Beckman-The-Zen-of-Expressing-State-The-State-Monad

We found Brian Beckman’s example overly complex, there’s quite a lot going on in that video (and the code needed de-obfuscating to read more like idiomatic C#), but it’s probably a good second or third example once the basics are understood. I wanted to look at a simpler example, and found this article: http://ertes.de/articles/monads.html#section-6 the explanation of the get and put methods show pretty clearly how the state monad works, and you should probably read and understand that section before carrying on here, as I’m going to jump straight on to show how the random number generator example looks in F# and C#.

C# implementation

First off, here’s an OO random number generator:

clip_image00125255B425255D

It is easy to add three random numbers together using an instance of the generator.

clip_image00225255B525255D

We can make NextShort purely functional, so that calling it has no side-effects, by passing the current state in, and returning the new state back out, along with the result:

clip_image00325255B525255D

To add three numbers, the state has to be kept around, and passed to each of the calls in turn. This obfuscates the main business logic.

clip_image00425255B425255D-1

Jumping straight to the punch line, using a state monad means that the boilerplate state management can be hidden away so that the algorithm reads as clearly as it did before:

clip_image00525255B425255D-1

 

This is a trivial example, but hopefully it makes the point. The rest of this blog post will show the implementation, along with the same in F#.

First off, we have a state class, which holds a Func that given a state, returns a tuple of the state and value. The F# example doesn’t bother with the State class, and I could have made the C# version simply use the Func everywhere, but it seemed clearer to me to have the wrapper object just to get a bit of context.

clip_image00625255B425255D-1

I created a non-generic State class as somewhere convenient to put the Get and Set methods:

clip_image00725255B425255D

The meat of the monad lives in extension methods in a StateEx class:

clip_image00825255B625255D

and SelectMany is implemented in terms of Bind the same as in Mike Hadlow’s example for the Identity and Maybe monads:

clip_image00925255B425255D

The computation to compose can be expressed in the GetRandom method:

clip_image01025255B425255D

This returns an operation (a Func in the State container) that invokes NextShort with the old state, and returns a computation with the new state and value.

The resulting composed operation using LINQ is described above, but this is what it looks like desugared:

clip_image01125255B425255D

F# implementation

The code for the random number generator and composed operation look similar to the C#:

clip_image012_thumb25255B125255D

This has the same problem of explicitly passing the state around. I can ruin the surprise again and show how using monads hides the explicit state passing:

clip_image01325255B425255D

This code makes use of F#’s support for monads via Computation Expressions. The code for the state builder is from here: http://www.navision-blog.de/2009/10/23/using-monads-in-fsharp-part-i-the-state-monad/

clip_image014_thumb25255B125255D

The Haskell example also had the >> operator, but of course F# uses that as the forward composition operator J

Instead I’ve implemented:

clip_image015_thumb25255B125255D

getRandom looks pretty much identical to the Haskell example:

clip_image016_thumb25255B125255D

 

The composed version using the StateBuilder is shown above, but here’s the desugared version:

clip_image01725255B425255D

The F# version is much clearer to read than the C#, in part due to the type inference meaning that the code has less pointless fluff. Having functions live in the module and not as part of a class aids that too (can simply call setState instead of State.Set in the C#)

Following through the simpler example really helped my understanding, and hopefully will see examples in the future where it will be sensible to use the state monad in my coding.

Visualising Sound

This is the first in a series of blog posts about visualising sound, in a similar way to an oscilloscope. The geek in me thought it’d be a fun thing to do, as well investigate different technological approaches to implementing the app (interesting as it’s a real-time performant app).

clip_image001

An oscilloscope has a single screen, which refreshes on a given time period, displaying a number of traces. The user can control that time period, as well as the gain (the y axis scale).

The top graph is essentially the same view as an oscilloscope containing a single trace, and a fixed time period (further blog posts may investigate varying time periods).

The bottom graph is a sliding window with a longer time period – this is the advantage of implementing an oscilloscope in code, we can create charts that aren’t really feasible in a classic CRT oscilloscope.

This series of blog posts will investigate implementing this screen using F# idioms, as well as using the Reactive Extension for .NET (RX Framework), and TPL Dataflow.

There are further things that could be implemented in future blog posts which may be interesting to see how the varying approaches look like:

· Trigger, or Trigger and hold: Only refresh the oscilloscope view once a certain trigger level has been reached. This is interesting as may want to include a certain number of immediately prior to the point where the trigger was set.

· Log many traces.

· Spectrum analyser/FFT.

· Digital filtering.

· Comparing traces, or more complicated math channels.

· Heatmap.

· Adjustable time offset (delay) – useful on the oscilloscope view to centre a waveform on the screen, or for when comparing two or more channels output.

F# Implementation

This blog post covers an F# implementation of the microphone, using asynchronous workflows. The graphing is being done by a nightly build download of Dynamic Data Display (D3). I’m really impressed with the performance, but it is quite difficult to customise.

The low-latency microphone code is from this article on gamedev (http://www.gamedev.net/topic/586706-slimdx-directsound-capture-microphone-stream-with-low-latency-example/).

The implementation of this is pretty simple; I’ll start from the detail out.

The inner loop of the program is an asynchronous workflow that reads from the buffer and returns a sequence:

clip_image003

Note the slightly-strange looking code:

let task = System.Threading.Tasks.Task<int>.Factory.StartNew(fun () -> WaitHandle.WaitAny(waitHandles))

F# has an Async.AwaitWaitHandle() method, which unfortunately only waits on a single handle. We want to wait on both handles so that we get notified when the buffer is full every 4096 instead of every 8192 bytes. With 4 bytes per sample, and a sample rate off 44KHz, this is equivalent to getting notified at an approximate rate of 40 times per second instead of 20 times per second.

I could have implemented Async.AwaitAnyWaitHandle() taking an array of WaitHandles, but looking at the code in the F# PowerPack, the code was quite complex. So, the code instead creates a new future to do the waiting and let us know which WaitHandle was set (this does mean that we’ve got the minor overhead of scheduling a new task to run on the task pool).

The Async.StartImmediate method ensures that the ProcessStream call is marshalled back onto the UI thread. It may be worth in the future looking at doing more of the data processing on a dedicated thread, leaving the UI thread free for drawing and user input.

The convertByteArrayToSequence is simple, it just iterates over the buffer in 4 byte chunks, and converts the values to floats, which it yields in the sequence:

clip_image004

The ProcessStream method looks like this:

clip_image006

For completeness, this is the Seq.sample module extension:

clip_image007

The nastiest bit of the code in ProcessStream is probably towards the end, where the windowedUnderlyingData is manipulated to ensure that the window only contains 1000 samples. It would be nice to do this in a non-imperative way, using Seq.windowed, but the problem is, is that the sequence we’re operating on is only the result of one buffer read operation, whereas windows etc. should operate over the whole data stream, and the sequences can’t be combined into a single sequence using yield! as they’re all generated asynchronously. Similarly, the buffer takes non-overlapping windows over the input sequence, and without taking a buffer over the whole stream, it may miss samples off the end of the sequence portions. Tomas Petricek has written a library, AsyncSeq, http://tomasp.net/blog/async-sequences.aspx which I may investigate in a later post.

The alternative to this would be to implement mailbox processors, having different agents for the buffering, windowing and sampling operations. I did start investigating this, but didn’t feel happy with them being pull-based (they don’t return data unless asked). I could have set them up to act more like a dataflow network, but it does seem to go against their intended use of managing state between different threads. I may revisit this in a future blog post.

I feel that even though F#’s does have nice features which helped to quickly implement the app, RX would probably be a better fit. I guess I won’t know until I implement it in RX and compare the differences.

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.