This project has moved and is read-only. For the latest updates, please go here.

F# Interface and Syntax

Mar 28, 2011 at 11:35 PM

Hello,

I am currently running numerical simulations and have been using F# for the task. F# has a flexible syntax that allows for expressive mathematical notation as well as having decent performance (with better performance possible via dropping into C++/CLI as needed). I like the clean OO API being used in Numerics, and I think that definitely is the way to go for the underlying numerical library. However, I was wondering if there would be any interest in creating a more "script-like" environment for the F# bindings so as to allow clean mathematical notation? A good example of this approach is the SciPy/NumPy/Matplotlib stack, where the core algorithms are in an OO hierarchy but they also expose a scripting environment called pylab which lets you write code in a quick and clean style similar to MATLAB/Octave.

There are really two parts to such an environment: more overloaded operators on the vector/matrix types and a library of functions which follows MATLAB-like naming conventions and allow easy generation of numerical code without getting into the depths of the OO hierarchy. For example, it'd be nice to be able to write a=zeros(3,3) and have a 3x3 matrix of zeros, without figuring out what storage type you want (as mentioned in a previous thread, 95% of the time you probably want Double.DenseMatrix or Complex.DenseMatrix). To that end I wanted to suggest reforming the F# bindings currently in the project to take advantage of the F# language's flexibility in creating a scripting environment.

Here is what some example code would look like:

 

let a = zeros(3,3)       // 3x3 zero matrix
let b = ones(3,3)        // 2x2 matrix of ones
let c = b.^3             // c=b*b*b
let d = b.^(-1)          // d=inverse(b) (should probably be pseudo-inv too?)
let e = b.*b             // e is all ones, elementwise multiplication
let f = b*b              // f is all threes, matrix multiplication
let g = b*b*2.0+b./2.0   // composition works fine 

 

The choice of operators can exactly copy MATLAB, since periods are usable as operators in F#. Here is a toy example of what a Kalman Filter would look like using such a binding:

 

let KalmanFilter (model:KFModel) (estimate:KFEstimate) (z:Matrix)=        
    let (x,P)=(estimate.x,estimate.P)
    let (H,R,Phi,Q)=(model.H,model.R,model.Phi,model.Q)

    // First propagate to the current time
    let x = Phi*x
    let P = Phi*P*Phi.T+Q

    // Now update using measurement
    let K = P*H.T*(H*P*H.T+R).^(-1)
    let res = z-H*x
    let x = x+K*res
    let P = P - K*H*P

 

The code is clean and readable by a mathematician/engineer without parsing parentheses. Of course, talk is cheap, so here is a partial implementation of a Matrix<T> wrapper type in F# which would allow the above code to run. I believe there is no easy way to do this without wrapping the C# class in a F# type unless the C# implementation of Matrix<T> is modified, as type extensions do not allow operator overloading and static type parameters aren't flexible enough. Of course, I'd love to hear a better way! Hopefully this is the right place for this, I wasn't sure what the codeplex conventions were for code pasting:

 

type Matrix<'T when 'T : (new:unit -> 'T) and 'T:struct and 'T :>System.ValueType and 'T :> System.IEquatable<'T> and 'T :> System.IFormattable> (inner:Generic.Matrix<'T>) =

    member this.inner with get() = inner

    member this.Item 
        with get(n,m)=inner.[n,m]
        and set(n,m) v=inner.[n,m]<-v
        
    override this.ToString()=
        inner.ToString()

    // Redefine transpose for brevity
    member this.T with get()=new Matrix<'T>(inner.Transpose())

    // Hermitian transpose
    member this.H with get()=new Matrix<'T>(inner.ConjugateTranspose())

    // Elementwise operations (matching MATLAB syntax)
    static member (.*) (a:Matrix<'T>,b:Matrix<'T>) = new Matrix<'T>(a.inner.PointwiseMultiply(b.inner))
    static member (./) (a:Matrix<'T>,b:Matrix<'T>) = new Matrix<'T>(a.inner.PointwiseDivide(b.inner))

    // Elementwise work on scalars (following MATLAB convention)
    static member (.*) (a:Matrix<'T>,b:'T) = new Matrix<'T>(a.inner.Multiply(b))
    static member (./) (a:Matrix<'T>,b:'T) = new Matrix<'T>(a.inner.Divide(b))
    

    /// Matrix scalar exponentiation for integer powers > 0, inverses
    static member (.^) (a:Matrix<'T>,b:int) = 
        if b= -1  then
            new Matrix<'T>(a.inner.Inverse())
        elif b<1 then
            raise (new NotImplementedException())
        else
            let mutable out = a.inner.Clone()
            for i in 0..b-2 do
                out<-out*a.inner
            new Matrix<'T>(out)

    /// Matrix-valued exponentiation (i.e. e^(matrix))
    static member (.^) (a:double,b:Matrix<'T>) =
        raise (new NotImplementedException()) // todo

    /// Matrix "right" division (i.e. conventional division)
    static member (/) (a:Matrix<'T>,b:Matrix<'T>) = 
        if a.inner.ColumnCount=a.inner.RowCount && b.inner.ColumnCount=b.inner.RowCount && b.inner.Rank()=b.inner.ColumnCount then
            new Matrix<'T>(a.inner*b.inner.Inverse())
        else
            raise (new NotImplementedException()) // todo: Should use pseudoinverse

    /// Matrix "left" division (i.e. mathemtically this performs a\b, but "\" is not a valid F# operator)
    static member ldiv (a:Matrix<'T>,b:Matrix<'T>) = 
        if a.inner.ColumnCount=a.inner.RowCount && b.inner.ColumnCount=b.inner.RowCount && a.inner.Rank()=a.inner.ColumnCount then
            new Matrix<'T>(a.inner.Inverse()*b.inner)
        else
            raise (new NotImplementedException()) // todo: Should use pseudoinverse    

    // Conventional operations. These are already overloaded in C# for the inner class
    static member (*) (a:Matrix<'T>,b:Matrix<'T>) = new Matrix<'T>(a.inner.Multiply(b.inner))
    static member (*) (a:Matrix<'T>,b:'T) = new Matrix<'T>(a.inner.Multiply(b))
    static member (*) (a:'T,b:Matrix<'T>) = new Matrix<'T>(b.inner.Multiply(a))
    static member (-) (a:Matrix<'T>,b:Matrix<'T>) = new Matrix<'T>(a.inner.Subtract(b.inner))
    static member (~-) (a:Matrix<'T>) = new Matrix<'T>(a.inner.Negate())
    static member (+) (a:Matrix<'T>,b:Matrix<'T>) = new Matrix<'T>(a.inner.Add(b.inner))
    static member (~+) (a:Matrix<'T>) = a

 

We'd also need the supporting library functions we used:

 

let zeros(n,m)=
    new Matrix(new Double.DenseMatrix(n,m))

let ones(n,m)=
    new Matrix(new Double.DenseMatrix(n,m,1.0))

let linspace(a,b,length)=
    let mutable out = zeros(1,length)
    for i in 0..length-1 do
        out.[0,i]<-a+double(i)*(b-a)/double(length-1)
    out

 

If this is something that interests anyone, I could start writing a more complete wrapper that wraps all the exposed members of Matrix<T>. Also, the current F# bindings (i.e. ofList, etc.) would need to be modified to support the F# wrapper class, as they are quite useful and idiomatic F#. Let me know what everyone thinks, and thanks for actually reading my giant wall of text/code! :)

P.S. while typing up this example I think I may have run into a bug. Apparently the code:

let g = MathNet.Numerics.LinearAlgebra.Double.DenseMatrix(2,2,3.0).Divide(2.0)
returns a zero matrix. Is that behavior intended?

Mar 29, 2011 at 5:04 PM

I would like to second this idea.  I'm not sure about implementation or anything, but F# seems ideally suited for numerics, so an F# interface would be great.  I can contribute cycles if needed.

Mar 31, 2011 at 3:30 PM
Edited Mar 31, 2011 at 3:47 PM
>let g = MathNet.Numerics.LinearAlgebra.Double.DenseMatrix(2,2,3.0).Divide(2.0)
>returns a zero matrix. Is that behavior intended?
I'll look in to it. Thanks for pointing it out. 
As for your suggestion, I like the idea but don't know enough about F# to really comment. Jurgen, what do you think?
Mar 31, 2011 at 3:48 PM
This discussion has been copied to a work item. Click here to go to the work item and continue the discussion.
Apr 5, 2011 at 12:19 PM

>let g = MathNet.Numerics.LinearAlgebra.Double.DenseMatrix(2,2,3.0).Divide(2.0)
>
returns a zero matrix. Is that behavior intended?

I wasn't able to duplicate the problem. Please try the most recent change set.

Apr 6, 2011 at 7:17 AM

Yeah my bad, all is well if I compile from source (was using beta 1).

Apr 15, 2011 at 7:36 PM

I've had some time to improve the F# interface I posted earlier:

https://gist.github.com/922192

There are two types defined here:
- FSMatrix is a wrapper around Matrix<T> that has the usual overloaded operators
- NumericsPlot is a class which exposes a simple plotting interface. You can pass it X and Y coordinates or just Y coordinates (X will be 1,2,3...). It accepts a list of points, array of points, Matrix<T> of points, or FSMatrix of points. Only works for integers/floats. When plotting matrices it plots each column as a separate data set. Note that plots are run in their own threads with their own event loops, so your program will not halt until the windows are closed (therefore if you want to just plot some stuff, you can exit after calling plot() and the program will stay open).

Here is an example of plotting a matrix/simple list:

open System
open FSNumerics
open MathNet.Numerics.LinearAlgebra

// Plot a list of points
let p = new NumericsPlot()
p.plot([1;2;3;4;5;4;3;2;1])

// Do an experiment, fill a matrix, print the results
let numRuns = 5
let numSamples = 50000
let data = zeros(numRuns,numSamples)

for i in 0..data.inner.RowCount-1 do
    for j in 0..data.inner.ColumnCount-1 do
        data.[i,j]<-sin(float(j)/10000.0-float(i*2))

p.plot(data.T) // Plots each column in the matrix in a separate plot, total of 5 plots

The resulting plots look like this:

http://i.imgur.com/hyeHp.png

May 27, 2011 at 10:55 PM

Hi kyon,

I like the succinct syntax of this new type. It's too bad we can use F# type augmentations to augment the existing types with the operators you've defined. I think that would be a neater solution rather than defining a whole new type. I'm not entirely sure we should add yet another specific F# type to the library; let me get in touch with some other people from the F# community and have them weigh in ...

(Sorry for the massive delay)

Jun 9, 2011 at 2:43 AM

Hi!
Jurgen asked me for some thoughts about the suggestion, so I played a bit with using Math.NET from F#. The wrapper definitely makes some common tasks nicer and there isn't really a way to get exactly the same extensions without using a wrapper. On the other hand, wrappers are (in principle) a bit painful. If you want to wrap everything, then it is a lot of code that must be maintained and if it is not complete, then people will need to explicitly wrap/unwrap the matrix themselves to get all the features (for example, the extension methods for Svd, Evd, GramSchmidt, etc. would have to be duplicated too).

You can get some of the features just by adding extension members to the type in the F# library. For example, here is a definition of the T and H members:

type MathNet.Numerics.LinearAlgebra.Generic.
  Matrix<'T when 'T : struct and 'T : (new : unit -> 'T) and 'T :> System.ValueType and 
                 'T :> System.IEquatable<'T> and 'T :> System.IFormattable> with
  // Redefine transpose for brevity
  member this.T = this.Transpose()

  // Hermitian transpose
  member this.H = this.ConjugateTranspose()

The additional operators are more tricky - they can be defined using "let" outside of the type, but then they cannot be overloaded (which means you'd need different operator for point-wise multiplication and multiplication by a scalar). One way to implement those is to use static member constraints, which essentially just says - if the type has "PointwiseMultiply" method, then it can use this operator:

let inline (.*) (a:^T) (b:^T) = 
  (^T : (member PointwiseMultiply : ^T -> ^R) (a, b))
let inline (./) (a:^T) (b:^T) = 
  (^T : (member PointwiseDivide : ^T -> ^R) (a, b))

let inline (.*|) (a:^T) (b:^S) = 
  (^T : (member Multiply : ^S -> ^R) (a, b))
let inline (./|) (a:^T) (b:^S) = 
  (^T : (member Divide : ^S -> ^R) (a, b))

The nice thing about this trick is that it will work for both matrices and vectors. You could even use the same operator for different types of operations (e.g. scalar and point-wise multiplication), but then the Matrix type would have to provide these two operations as some overloaded method with the same name (it could be something like "FSharpOverloadedMultiply", but that would look ugly from C#). It is not as nice as what the wrapper can provide, but it looks decent:

let zeros n m = new Double.DenseMatrix(n,m)
let ones n m = new Double.DenseMatrix(n,m,1.0)

let v = vector [ 1.0; 2.0 ]
v .* (v .*| 2.0)

let m = matrix [ [ 1.0; 2.0]; [ 3.0; 4.0 ] ]
m.H .* ((ones 2 2) .*| 2.0)

Regarding charting - there is an F# sample called FSharpChart that implements a pretty complete F# wrapper for the .NET WinForms charting libraries. The code is available as a sample with permissive license, so it should be fine to include it in the F# library for Math.NET Numerics (or in some other part of Math.NET). That would be actually quite nice, because currently it is "just" a sample.

I don't have any clear suggestion on what is better design - adding F# specific functionality using extensions and additional operators doesn't give you all the features that a wrapper can give, but wrappers have quite a few issues. I guess it is also perfectly fine to have the wrapper in some additional script file that people can reference if they want, but maybe not in the standard library.

Hope this helps!
Tomas