## Tuesday, 27 July 2010

### F#unctional Londoners meetup lecture (28th July 2010)

Zach Bray of Trayport and Jon Harrop of Flying Frog Consultancy Ltd. will be presenting lectures at Skills Matter eXchange London (UK) at 6:30pm on Wednesday 28th July 2010.

Many thanks to Carolyn Miller and Phil Trelford for organizing F#unctional Londoners Meetup Group, an excellent series of events!

## Saturday, 17 July 2010

### Histogram equalization

The F#.NET Journal just published an article about image processing and the new charting functionality in .NET 4:

"Over- and under-exposed images have their intensity distributions skewed to the high or low end of the range. Histogram equilization is one technique to combat this effect digitally by spreading out the distribution of intensities in an image via the intensity histogram. This article describes a program for image enhancement using histogram equalization with an interactive WPF-based GUI using the new charting functionality in .NET 4 to visualize the intensity histograms in real time..."

## Monday, 12 July 2010

### F# vs Mathematica: an even faster pricer for American options

The previous post described a performance-critical pricer for American options. Upon reflection, the optimized Mathematica code had introduced many redundant computations in order to vectorize the code because vectorization offers considerable performance improvements in such a slow language. Our simple translation of the optimized Mathematica had not taken advantage of this by removing these redundant computations.

This new optimized and parallelized F# solution is now a whopping 960× faster than the original Mathematica code from Sal Mangano's Mathematica Cookbook:

```let americanPut kk r sigma tt =
let sqr x = x * x
let ikk, a, nn, mm, tt0 = 1.0 / kk, 5.0, 1000, 200, sqr sigma * tt / 2.0
let k, h, s = 2.0 * r / sqr sigma, 2.0 * a / float nn, tt0 / float mm
let alpha, k0, k1 = s / sqr h, 0.5 * (k - 1.0), 0.25 * sqr (k + 1.0)
let x i = float i * h - a
let ex = Array.init (nn+1) (fun i -> exp(x i))
let ek0x = Array.init (nn+1) (fun i -> exp(k0 * x i))
let pp0 i = max 0.0 (kk - kk * ex.[i])
let ek0a, k1s, a' = exp(k0 * a), k1 * s, 2.0 * alpha - 1.0
let mutable u = Array.init (nn+1) (fun i -> ek0x.[i] * pp0 i * ikk)
let mutable u' = Array.create (nn+1) 0.0
for j=0 to mm-1 do
let ek1sj = exp(k1s * float j)
let k2 = alpha * ikk * ek0a * ek1sj
u'.[0] <- alpha * u.[1] - a' * u.[0] + k2 |> max (ek0x.[0] * ek1sj * (1.0 - ex.[0]))
let inline get (u: _ []) i m = max m (alpha * (u.[i+1] + u.[i-1]) - a' * u.[i])
for i=1 to nn/2 do
u'.[i] <- get u i (ek0x.[i] * ek1sj * (1.0 - ex.[i]))
for i=nn/2 to nn-1 do
u'.[i] <- get u i 0.0
let u'' = u
u <- u'
u' <- u''
let k3 = kk * exp(-k1 * tt0)
Seq.mapi (fun i u -> kk * ex.[i], k3 / ek0x.[i] * u) u

Array.Parallel.init 1000 (fun i ->
americanPut (float(i+1)) 0.05 0.14 1.0)```

Note that the hardcoded constants nn and mm have been increased, sigma has been decreased to keep alpha just below 0.5 and strikes from 1 to 1,000 are computed in order to make the F# code take a significant fraction of the second (when the Mathematica takes several minutes!).

## Friday, 9 July 2010

### F# vs Mathematica: fast pricer for American options

UPDATE: A new F# solution that is 960× faster than the original Mathematica is available here.

Another example from Sal Mangano's Mathematica Cookbook, this time taken from Andreas Lauschke's example, is a "fast" pricer for American options. The relevant section of the book stresses the importance of performance in this context and the Mathematica code presented was heavily optimized by experts. Here is their code:

Specifically, this function has been optimized by pushing computational burden from the general-purpose term rewriting language of Mathematica onto the optimized C routines in its standard library and then compiling the high-level code into lower-level bytecode that is interpreted more efficiently, giving another 5× speedup.

However, this general approach to the optimization of Mathematica code has the unfortunate side effect of "foresting": introducing unnecessary intermediate data structures because built-in operations over them are faster than writing a direct solution in such a slow language. Consequently, a direct solution in a compiled language will usually be orders of magnitude faster. In this case, we found that the following translation to F# is 64× faster than Sal's original Mathematica benchmark:

```let americanPut kk r sigma tt =
let sqr x = x * x
let a, nn, mm, tt0 = 5.0, 100.0, 20, sqr sigma * tt / 2.0
let k, h, s = 2.0 * r / sqr sigma, 2.0 * a / nn, tt0 / float mm
let alpha, xn = s / sqr h, int(2.0 * a/h + 0.5) + 1
let x i = float i * h - a
let ss = Array.init xn (fun i -> kk * exp(x i))
let k0, k1 = 0.5 * (k - 1.0), 0.25 * sqr (k + 1.0)
let k1s = k1 * s
let pp0 i = max 0.0 (kk - ss.[i])
let rec run (u: _ []) (u': _ []) j =
if j=mm then u else
let f i =
let xi = -a + float i * h
if xi >= 0.0 then 0.0 else
exp(k0 * xi + k1s * float j) * (1.0 - exp xi)
u'.[0] <- alpha * u.[1] - (2.0 * alpha - 1.0) * u.[0] +
alpha / kk * exp(k0 * a + k1s * float j) |> max (f 0)
for i=1 to u.Length-2 do
u'.[i] <- alpha * (u.[i+1] + u.[i-1]) -
(2.0 * alpha - 1.0) * u.[i] |> max (f i)
run u' u (j+1)
let u i = exp(k0 * x i) * pp0 i / kk
let u = run (Array.init xn u) (Array.create xn 0.0) 0
ss, Seq.mapi (fun i u -> kk * exp(-k1 * tt0 - k0 * x i) * u) u

Array.Parallel.init 10000 (fun strike ->
let ss, ps = americanPut (float(strike+1)) 0.05 0.4 1.0
Seq.zip (Seq.take 60 ss) (Seq.take 60 ps))```

The superior performance of the F# solution comes from two main benefits:

• The F# program is compiled all the way down to machine code before being executed whereas the Mathematica code is interpreted. The intermediate data structures are then replaced by simple function calls. This alone makes the F# solution over 10× faster than the original Mathematica.

• F# inherits the benefits of a highly-optimized multicore-capable run-time from the .NET platform. This allows us to use fine-grained parallelism to improve performance even further for another 6.1× speedup on this 8-core machine.

## Thursday, 8 July 2010

### F# vs Mathematica: parametric plots

Another example from Sal Mangano's Mathematica Cookbook (p.520) is an elegant little Mathematica program that uses a crude numerical integrator to plot the trajectory of a differential equation representing the populations of predators (foxes) and prey (rabbits):

This example, which was derived from this Wolfram Demonstrations sample, really plays to the strengths of the Mathematica language, not least its enormous standard library of graphing libraries.

However, .NET 4 does include new charting functionality so it is interesting to see how elegantly this can be written in F#. The following F# program implements the same functionality when evaluated interactively:

```#r "System.Windows.Forms.dll"
#r "System.Windows.Forms.DataVisualization.dll"
#r "WindowsBase.dll"
#r "WindowsFormsIntegration.dll"
#r "PresentationCore.dll"
#r "PresentationFramework.dll"
#r "System.Xaml.dll"

open System.Windows.Forms
open System.Windows.Forms.DataVisualization.Charting

let trajectory g k t =
let evolve(r, f) =
let dtrf = 0.0001 * r * f
r + (1.0 - r/k)*r*g - dtrf, dtrf + (1.0 - g)*f
Seq.scan (fun s _ -> evolve s) (50.0, 10.0) {1..t}

use series = new Series(ChartType=SeriesChartType.Line)
for x, y in trajectory 0.02 5e2 1500 do
use area = new ChartArea()
area.AxisX.Title <- "Rabbits"
area.AxisX.Minimum <- 0.0
area.AxisY.Title <- "Foxes"
use chart = new Chart(Dock=DockStyle.Fill)
use form = new Form()
form.Show()```

When evaluated, this F# program produces the following output:

This program will be the basis for a future F#.NET Journal article that introduces user interaction. Although Mathematica also makes user interaction easy as well, its relatively poor performance often makes interactive programs sluggish whereas the relevant code in our F# version of this program runs 64× faster than the original Mathematica and, consequently, produces a fluid user interface.

## Wednesday, 7 July 2010

### F# vs Mathematica: red-black trees

Mathematica is an expensive commercial application sold as software for "interactive technical computing". In particular, the product centers around the Mathematica programming language which is a very powerful and dynamic language based upon term rewriting with excellent graphical capabilities and a wealth of useful functionality in its standard library. However, Microsoft's new F# programming language not only provides most of the useful functionality found in Mathematica but has many other benefits, not least that it is completely free!

We have previously compared F# and Mathematica in the context of 2D vector graphics. This post is the first in a series comparing these languages. Specifically, this post compares data structures. The F# programming language is designed for efficient compilation to native code whereas the Mathematica language is generally evaluated via term rewriting which is, essentially, a very inefficient form of interpretation. Consequently, F# code is typically 10-1,000× faster than Mathematica code unless computational burden can be shifted onto functions written in other languages (usually those in the vast standard library).

The latest Mathematica-related book, Mathematica Cookbook by Sal Mangano, contains many examples including a translation of Chris Okasaki's immutable red-black trees from Haskell to Mathematica augmented with a function to remove an element from a tree. Unfortunately, the new function written by Sal Mangano is incorrect and silently corrupts the shape of the tree. Thus, we shall focus on the subset of the code presented in Sal's book that is correct. This code may be written as follows in F#:

```type color = Red | Black
type 'a t = Node of color * 'a * 'a t * 'a t | Leaf

let balance = function
| Black, z, Node (Red, y, Node (Red, x, a, b), c), d
| Black, z, Node (Red, x, a, Node (Red, y, b, c)), d
| Black, x, a, Node (Red, z, Node (Red, y, b, c), d)
| Black, x, a, Node (Red, y, b, Node (Red, z, c, d)) ->
Node (Red, y, Node (Black, x, a, b), Node (Black, z, c, d))
| x -> Node x

let insert x s =
let rec ins = function
| Leaf -> Node (Red, x, Leaf, Leaf)
| Node (color, y, a, b) as s ->
if x < y then balance (color, y, ins a, b)
elif x > y then balance (color, y, a, ins b) else s
match ins s with
| Node (Red, y, a, b) -> Node (Black, y, a, b)
| t -> t```

The Mathematica code is significantly more complicated (714 vs 1514 bytes) for two main reasons:

• Mathematica has no types so it cannot express a type with an ordering, so a new rbTree data structure is created with an ordering and that must then be boxed and unboxed by hand.
• Mathematica prohibits top-level or-patterns so the same pattern is repeated four times in the balance function.

The simpler F# solution is also 20× faster at inserting 100k elements and 100× faster at computing the range of depths in the resulting tree.

In the next post, we'll start to compare the graphical capabilities of the two languages.

## Monday, 5 July 2010

### Happy numbers

Consider the sequence obtained by summing the squares of the digits of a number and then repeating with the result. If this sequence ends with the number one then the numbers in the sequence are Happy Numbers.

One of the challenges on Rosetta code is to compute the first eight Happy Numbers. The current F# solution weighs in at 26 lines of code. Let's see if we can do better...

We begin by writing a next function that computes the next happy number after the given number:

```> let rec next (m: Set<_>) n i =
if i=1 then n else
if m.Contains i then next Set.empty (n+1) (n+1) else
Seq.sumBy (fun d -> pown (int d - int '0') 2) (string i)
val next : Set<int> -> int -> int -> int```

Note the use of the string function to convert an int into its string representation and then the use of the int function to convert each char back into an int.

Finally, we use the List.scan function to build a list of accumulated results as the next function finds each of the happy numbers in turn:

```> List.scan (fun n _ -> next Set.empty (n+1) (n+1)) 1 [1..7];;
val it : int list = [1; 7; 10; 13; 19; 23; 28; 31]```

The scan function is relatively new and, consequently, not as widely appreciated as perhaps it should be.

Our whole program weighs in at only 6 lines of code, much shorter than the original!

## Sunday, 4 July 2010

### Lorenz attractor

The Lorenz attractor is a fractal derived from the trajectory of a 3-dimensional dynamical system that exhibits chaotic flow. This blog post describes a 35-line program that computes trajectories of this attractor and visualizes them as a beautiful whisp using Windows Presentation Foundation:

The program is as follows:

```> #r "PresentationCore.dll";;
> #r "PresentationFramework.dll";;
> #r "WindowsBase.dll";;
> #r "System.Xaml.dll";;

> let trajectory f (x, y, z) dt n =
let s, b, p = 10.0, 8.0 / 3.0, 28.0
let mutable x = x
let mutable y = y
let mutable z = z
let a = Array.zeroCreate n
for i=0 to n-1 do
a.[i] <- f (x, y, z)
x <- x + s * (y - x) * dt
y <- y + (x * (p - z) - y) * dt
z <- z + (x * y - b * z) * dt
a;;
val trajectory :
(float * float * float -> 'a) ->
float * float * float -> float -> int -> 'a []

> open System.Windows;;

> let whisp() =
let rand = System.Random()
let f x = x + pown (rand.NextDouble() - 0.5) 2
let x, y, z = f 10.0, f 0.0, f 20.0
let xys =
let f(x, y, z) = Point(20.0 + x, 25.0 + y - z + 50.0)
trajectory f (x, y, z) 0.003 2000
let line_to (xy: Point) =
(Media.LineSegment(xy, true) :> Media.PathSegment)
Media.PathGeometry[ Media.PathFigure(xys.[0], Seq.map line_to xys, false) ];;
val whisp : unit -> Media.PathGeometry

> do
let group = Media.GeometryGroup()
for i=1 to 100 do
let brush = Media.SolidColorBrush Media.Colors.Red
let path = Shapes.Path(Data=group, Stroke=brush, StrokeThickness=0.005)
let box = Controls.Viewbox(Child=path, Stretch=Media.Stretch.Uniform)
let window = Window(Content=box, Title="Lorenz attractor")
(Application()).Run window
|> ignore;;```

The trajectory function generates an array derived from the coordinates along a trajectory with the given starting position. The whisp function computes a trajectory with a slightly randomized starting position and generates a WPF PathGeometry. Finally, the main program generates 100 whisps and visualizes them.

### Rigid body dynamics demo

The F#.NET Journal just made another downloadable demo freely available. This demo is from the article Rigid body dynamics (15th January 2010).

The article walks through the design and implementation of a 250-line program that simulates the dynamics of a collection of solid balls as they bound around a scene, using Windows Presentation Foundation for visualization!

## Friday, 2 July 2010

### Diff

The F#.NET Journal just published an article about text processing:

"The unix program diff identifies differences between text files, line by line. This tool is most useful for comparing two versions of a program. This article develops a simple implementation of the diff tool based upon the longest common subsequence (LCS) algorithm with a graphical user interface written using Windows Presentation Foundation to visualize the difference between two text files..."