Tutorial 10: Generic programming using quotations

The code for this tutorial comes in several files.
First: GenericPhysics.fs

#light

open System
open Microsoft.FSharp.Linq.QuotationEvaluation

#nowarn "57"


type ps<'Vec> = class
val pos:'Vec
val speed:'Vec
new(pos, speed) = { pos = pos; speed = speed }
end


type sa<'Vec> = class
val speed:'Vec
val accel:'Vec
new(speed, accel) = { speed = speed; accel = accel }
end


let mk_gravity scale_func (up: 'Vec): 'Vec =
let q =
<@
let (*) = %scale_func in -9.81 * up
@>
q.Eval()


let mk_spring scale_func: float -> 'Vec -> 'Vec =
let q =
<@
let (*) = %scale_func in fun k v -> -k * v
@>
q.Eval()


let mk_drag scale_func: float -> 'Vec -> 'Vec =
let q =
<@
let (*) = %scale_func in fun k v -> -k * v
@>
q.Eval()


let mk_euler_implicit plus_func scale_func =
let q =
<@
let (+) = %plus_func in
let (*) = %scale_func in
fun accel_func pos speed delta ->
let speed2 = speed + delta * (accel_func pos speed)
let pos2 = pos + delta * speed2
ps<_>(pos2, speed2)
@>
q.Eval()


let inline adapt (f: 'Vec -> 'Vec -> 'Vec) =
fun pos speed ->
let accel = f pos speed in
sa<_>(speed, accel)


let mk_rk4_q plus_func scale_func =
<@
fun deriv_func pos speed delta ->
let pf = %plus_func
let sf = %scale_func

let extrapolate x delta dx =
pf x (sf delta dx)

let half_delta = 0.5 * delta

let sa1: sa<_> = deriv_func pos speed
let sa2 = deriv_func (extrapolate pos half_delta sa1.speed) (extrapolate speed half_delta sa1.accel)
let sa3 = deriv_func (extrapolate pos half_delta sa2.speed) (extrapolate speed half_delta sa2.accel)
let sa4 = deriv_func (extrapolate pos delta sa3.speed) (extrapolate speed delta sa3.accel)

let sixth_delta = delta / 6.0

let update old v1 v2 v3 v4 =
pf old (sf sixth_delta
(pf v1
(pf v2
(pf v2
(pf v3
(pf v3
v4))))))

let pos1 = update pos sa1.speed sa2.speed sa3.speed sa4.speed
let speed1 = update speed sa1.accel sa2.accel sa3.accel sa4.accel

ps<'Vec>(pos1, speed1)
@>


let mk_rk4 plus_func scale_func =
let q = mk_rk4_q plus_func scale_func
q.Eval()


let simulate intg_func (t0: float) t_fin delta pos0 speed0 =
let rec repeat t0 pos0 speed0 =
if t0 > t_fin then (pos0, speed0)
else
let t1 = t0 + delta
let ps: ps<_> = intg_func pos0 speed0 delta
repeat t1 ps.pos ps.speed
repeat t0 pos0 speed0

Second file: FloatPhysics.fs

#light

let gravity = GenericPhysics.mk_gravity <@ fun x y -> x * y @> 1.0
let drag = GenericPhysics.mk_drag <@ fun x y -> x * y @>
let spring = GenericPhysics.mk_spring <@ fun x y -> x * y @>
let rk4 = GenericPhysics.mk_rk4 <@ fun x y -> x + y @> <@ fun x y -> x * y @>
let euler = GenericPhysics.mk_euler_implicit <@ fun x (y:float) -> x + y @> <@ fun x y -> x * y @>

Third file: Util.fs

#light

open System

let run_time_func (func: unit -> 'a) =
let watch = new System.Diagnostics.Stopwatch()
watch.Start()
let res = func ()
watch.Stop()
let diff0 = float watch.ElapsedMilliseconds / 1000.0
(diff0, res)


let map_parallel func items =
let tasks =
seq {
for i in items -> async {
return (func i)
}
}
Async.Run (Async.Parallel tasks)

Last file: Main.fs

#light

open System
open FloatPhysics
open Util

let inline accel pos speed =
drag 1.0 speed + spring 100.0 pos + gravity

let ef = (FloatPhysics.euler accel)
let rk4f = (rk4 (GenericPhysics.adapt accel))

let run euler_func s0 =
let t0 = 0.0
let t_fin = 1000.0
let delta = 0.005

let run_timed map_func =
let n_runs = 2
let times =
[
for i in 1..n_runs ->
let (run_time, res) =
run_time_func (fun () ->
s0 |> map_func(fun (x, y) ->
GenericPhysics.simulate euler_func t0 t_fin delta x y))
run_time
]

let max = Seq.max times
let min = Seq.min times
let avg = (Seq.fold (fun x y -> x+y) 0.0 times) / float n_runs
(min, max, avg)

let single_min, single_max, single_avg = run_timed List.map
let multi_min, multi_max, multi_avg = run_timed map_parallel

printfn "Single (min, max, avg): %f %f %f" single_min single_max single_avg
printfn "Multi (min, max, avg): %f %f %f" multi_min multi_max multi_avg


let n_bodies = 2

let s0_float = [ for s in 1..n_bodies -> (float s, 0.0) ]
run ef s0_float
run rk4f s0_float

ignore(Console.ReadKey(true))

If you want to compile these sources using Visual Studio, you need to know that the order in which the files appear in the project is important. If a file B uses types or values from a file A, then A must be placed higher in the list of files of the project.

You will also need to add references to the following DLLs: FSharp.PowerPack and FSharp.PowerPack.Linq.

Let us now take a deeper look at the new features of F# used in the code.

open Microsoft.FSharp.Linq.QuotationEvaluation

#nowarn "57"

The main item of this Tutorial is quotations. These represent pieces of code which can be constructed at compile time, changed dynamically at run time, and recompiled and executed dynamically. This last part is still experimental, and requires the use of a specific module.
The "#nowarn" directives tells the compiler not to warn us about this feature being experimental.

let mk_gravity scale_func (up: 'Vec): 'Vec =
let q =
<@
let (*) = %scale_func in -9.81 * up
@>
q.Eval()

The interesting piece here is the part between <@ and @>. The entire thing is called a quotation. It must contain a valid segment of F# code. The F# compiler will parse the code in this quotation, and build an object representing this code. More precisely, the compiler will generate an abstract syntax tree (AST). Note that this object does not in itself constitute executable code.
You may be wondering what "%scale_func" means in this context. The "%" operator specifies that another AST should be inserted there. "scale_func" is the variable which should contain the AST to insert.
To clarify this, look at the call to mk_gravity in FloatPhysics.fs:

let gravity = GenericPhysics.mk_gravity <@ fun x y -> x * y @> 1.0

The resulting quotation would therefore look as follows:

<@
let (*) = fun x y -> x * y in -9.81 * up
@>

which is has the same effect as:

<@
-9.81 * up
@>

An AST generated from a quotation can be compiled using method "Compile", or compiled and executed using method "Eval".
Had I written "q.Eval()" instead, "mk_gravity" would have returned a value of type "unit -> 'Vec", i.e. a function taking no argument and returning a vector when called. We may just as well returned the vector directly, which we can do by using method "Eval" instead.

If we take a step back and look at what we have just achieved, we will see that we have constructed a rather complicated way of perform the product of a vector by a scalar.

Let us look at the next function for a more interesting example.

let mk_spring scale_func: float -> 'Vec -> 'Vec =
let q =
<@
let (*) = %scale_func in fun k v -> -k * v
@>
q.Eval()

While I wrote this code I repeatedly found myself wondering "What did I just write?". Happily, the compiler can answer that question.This function takes an AST which should represent a function taking a float and a vector (in curried form) and returning a vector. The "mk_spring" function will then return a function taking a float and a vector and returning a vector. In other words, it will "unquote" the quotation, which is simply done by evaluating the AST. Before doing so, it inserts the quoted form of the function performing the product of a vector by a scalar into q, which in the present case has the same effect as:

<@
let (*) = fun x y -> x * y in fun k v -> -k * v
@>

Or in simpler terms:

<@
-k * v
@>

We have found yet another way to perform the product of a vector (which is actually a float in this simple example) by a scalar (another float, but not to be mistaken with the vector). The "mk_spring" function differs from "mk_gravity" in the way that its parameters are expected to vary. Gravity is assumed to be constant, and so is the upward direction, denoted by the "up" parameter of "mk_gravity". The stiffness coefficient and the position of the body, on the other hand, change across calls.

Using quotations to achieve genericity can be seen as some kind of hybrid between the functional approach described in Tutorial 8 and the approach based on static constraints of Tutorial 9. It shares with Tutorial 8 its use of function parameters to specify the implementation of operations on vectors. It share with Tutorial 9 its way of emitting specific code for each concrete type of vector.

I had hoped this would result in more efficient code, but that is unfortunately not currently the case. The quotation evaluation functionality is still experimental, and it may be possible that future versions of F# might produce faster code.

Update
I've come across an entry in Jomo Fisher's blog about Linq which makes me wonder if my code is really right. Jomo has used Linq in a setting where performance is truly impressive. Why I am getting such bad results here? Further investigation needed...
End of Update

The point of the approach described here is not so much what it achieves, but what possibilities it opens. For instance, it is possible to implement aggressive optimization techniques that the F# compiler will probably never support, such as rearranging code to avoid creating objects for intermediate values, or turning pure computations into computations with local side effects. Obviously, it would not be realistic to support all features of F#, but a subset dedicated to computations would probably suffice. Using Quotations.Expr from F# and Reflection.Emit from .NET, it seems very feasible to write a specialized optimizing compiler for a subset of F#.
Other interesting potential uses of quotations for games include programs to run on the GPU (shaders), or on any other kind of peripheral, for that matter (game controllers, head-on displays...).