Computer Science Experimentation

Monday, September 27, 2010

F# parallel programming using Agents

Consider a parallel program as a series of independent components that send and receive messages.This is referred as the agent model.

When an actor receives a message, the message is placed in a queue until the agent is ready to process it. When an agent process a message, it makes a decision about what to do with it based on its internal state and contents of the message. It might do processing, send a reply, send a new message for a different agent, etc.

F# provides the MailboxProcessor class for agent processing using Asynchronous Programming or more specific, Asynchronous Workflows.

Asynchronous Programming describes programs and operations that once started are executed in background and terminate at some later time.

Asynchronous Workflow allow you to perform asynchronous operations without the need for explicit callbacks. So you can write the code as if it were using synchronous execution, but actually, the code will execute asynchronously, supending and resuming as asynchronous operation complete.

In Asynchronous Programming, you want to avoid blocking of threads. Threads are expensive because allocate a 1 MB stack and other operational system implications.
In performance critical code, it is important to keep the number of blocked threads low. Ideally, you will only have as many as you have processors and you will have no blocked threads.

There are several topics related to the usage of MailboxProcessor. In this document I will analyze the handling of agent's persisted variables.

The code for the first agent example is:

let Act_RandomAccum= 
    MailboxProcessor.Start(fun (i:MailboxProcessor) ->
        let rec loop state =
            async {
                    let! r= i.Receive()
                    let z=new Random()
                    let v=z.NextDouble()
                    printfn "state=%f" (state+v)
                    return! loop(state+v)}
        loop 0.0)
(*
state=0.761920
state=1.613361
state=1.717945
state=1.912050
...
*)

Let's create an function that will run continuously and it will post messages to agents every 3 seconds (asynchronously). You can cancel it by executing Async.CancelDefaultToken().

let Act_clk=
    async {
            while true do
                do Act_RandomAccum.Post(1)
                //do Act_3States.Post(1)
                //do Act_3States_2.Post(1)
                do! Async.Sleep(3*1000)
           }

let cancelHandler(ex : OperationCanceledException)=
    printfn "Task canceled"

Async.TryCancelled(Act_clk, cancelHandler) 
|> Async.Start

Notice how variables are keep by being a input to the loop.
This is not covinient when you want to preserve several states.

The next snippet uses a record to preserve several variables. Change the comment lines in Act_clk to send messages to the following agent.

type States= {mutable s1:float; mutable s2:float; mutable s3:float; mutable sum:float; mutable Status: string}

let States1={s1=0.0;s2=0.0;s3=0.0;sum=0.0;Status="good"}

let Act_Random3States= 
    MailboxProcessor.Start(fun (i:MailboxProcessor) ->
        let rec loop (s:States) =
            async {
                    let! r= i.Receive()
                    let z=new Random()
                    let v=z.NextDouble()
                    s.s3<-s.s2; s.s2<-s.s1
                    s.s1<-v
                    s.sum<-s.s1+s.s2+s.s3
                    printfn "states,sum=%f %f %f %f" s.s1 s.s2 s.s3 s.sum
                    return! loop(s)}
        loop States1)
(* 
states,sum=0.023424 0.000000 0.000000 0.023424
states,sum=0.046045 0.023424 0.000000 0.069468
states,sum=0.068666 0.046045 0.023424 0.138134
states,sum=0.568861 0.068666 0.046045 0.683572
states,sum=0.591482 0.568861 0.068666 1.229009
states,sum=0.614103 0.591482 0.568861 1.774447
states,sum=0.636724 0.614103 0.591482 1.842310
...
*)

The last example is more powerfull. The states are preserved in a class. Objects (as the random one) can be stored and don't need to be recreated in every execution.
Instance and static function can be defined (as fShift and fSum) to be used by the agent.
In other words, the agent defition is encapsulated in a class and its execution is done in the MailboxProcessor.

type StatesClass(s1:float,s2:float,s3:float,status:string)=
    let mutable _s1=s1 
    let mutable _s2=s2 
    let mutable _s3=s3 
    let mutable _status=status 
    let r=new Random()
    let mutable sum= StatesClass.fSum _s1 _s2 _s3  

    member this.Sum with get()=sum and set x=sum<-x
    member this.S1 with get()= _s1 and set x=_s1<-x
    member this.S2 with get()= _s2 and set x=_s2<-x
    member this.S3 with get()= _s3 and set x=_s3<-x
    member this.Status with get()= _status and set x=_status<-x
    member this.rndObj=r    
    static member fSum x y z=x+y+z
    member this.fShift()=
        this.S3<-this.S2; this.S2<-this.S1

let sc1=StatesClass(s1=0.0,s2=0.0,s3=0.0,status="ok")

let Act_3States_2= 
    MailboxProcessor.Start(fun (i:MailboxProcessor) ->
        let rec loop (s:StatesClass) =
            async {
                    let! r= i.Receive()
                    s.fShift()
                    let v=s.rndObj.NextDouble()
                    s.S1<-v
                    s.Sum <- StatesClass.fSum s.S1 s.S2 s.S3
                    printfn "states,sum=%f %f %f %f" s.S1 s.S2 s.S3 s.Sum
                    return! loop(s)}
        loop sc1)
(*
states,sum=0.396605 0.000000 0.000000 0.396605
states,sum=0.029249 0.396605 0.000000 0.425854
states,sum=0.197748 0.029249 0.396605 0.623602
states,sum=0.863481 0.197748 0.029249 1.090479
states,sum=0.641731 0.863481 0.197748 1.702961
states,sum=0.867453 0.641731 0.863481 2.372665
states,sum=0.880071 0.867453 0.641731 2.389254
...
*)

Thursday, September 23, 2010

Large-scale nonlinear optimization for F#

This blog shows how to use Ipopt (https://projects.coin-or.org/Ipopt) in F#.
Ipopt (Interior Point OPTimizer) is a software package for large-scale nonlinear optimization.
It is designed to find (local) solutions of mathematical optimization problems of the from
   min f(x)
x in R^n
s.t.    g_L <= g(x) <= g_U
        x_L <= x <= x_U

where f(x): R^n --> R is the objective function, and g(x): R^n --> R^m are the constraint functions. The vectors g_L and g_U denote the lower and upper bounds on the constraints, and the vectors x_L and x_U are the bounds on the variables x. The functions f(x) and g(x) can be nonlinear and nonconvex, but should be twice continuously differentiable. Note that equality constraints can be formulated in the above formulation by setting the corresponding components of g_L and g_U to the same value.

The binaries can be downloaded from http://www.coin-or.org/download/binary/Ipopt/. I am using the library Ipopt38.dll.
The interface for .NET can be downloaded from http://code.google.com/p/csipopt/.
Following the instruction on README.txt, you can create a class library Cureos.Numeric.dll by running the batch script build_library.bat.
Based on the C# example, provided by csipopt, I created the following F# code. I created a type, following the example, but it is not really necessary.

#I @"C:\Users\caxelrud\Documents\Visual Studio 2010\Projects\hs071_fsx"
#r "Cureos.Numerics.dll"
//Need Ipopt38.dll in path

open System
open System.IO
open Cureos.Numerics
open System.Runtime.InteropServices

type HS071=
    val public n : int
    val public m : int
    val public nele_jac : int
    val public nele_hess : int
    val public x_L : double[]
    val public x_U : double[]
    val public g_L : double[]
    val public g_U : double[]
    //set the number of variables and allocate space for the bounds
    //set the values for the variable bounds using n, x_L, x_U
    //Set the number of constraints and allocate space for the bounds using m, g_L, g_U
    //Set the number of nonzeros in the Jacobian of the constraints using nele_jac
    //Set the number of nonzeros in the Hessian of the Lagrangian (lower or upper triangual part only) using nele_hess

    new ()= {n=4;m=2;nele_jac = 8;nele_hess = 10;x_L =[1.0; 1.0; 1.0; 1.0 ];x_U = [5.0; 5.0; 5.0; 5.0];g_L=[25.0; 40.0];g_U = [Ipopt.PositiveInfinity; 40.0]}
 
    member public this.eval_f (n:int, x:double[], new_x:bool, [] obj_value:byref)=
        obj_value <- x.[0]*x.[3]*(x.[0]+x.[1]+x.[2])+x.[2]
        true

    member public this.eval_grad_f (n:int, x:double[], new_x:bool, grad_f:double[])=
        grad_f.[0] <- x.[0] * x.[3] + x.[3] * (x.[0] + x.[1] + x.[2])
        grad_f.[1] <- x.[0] * x.[3]
        grad_f.[2] <- x.[0] * x.[3] + 1.0
        grad_f.[3] <- x.[0] * (x.[0] + x.[1] + x.[2])
        true

    member public this.eval_g (n:int, x:double[], new_x:bool, m:int, g:double[])=

        g.[0] <- x.[0] * x.[1] * x.[2] * x.[3];
        g.[1] <- x.[0] * x.[0] + x.[1] * x.[1] + x.[2] * x.[2] + x.[3] * x.[3];
        true;

    member public this.eval_jac_g (n:int, x:double[], new_x:bool, m:int, nele_jac:int, iRow:int[], jCol:int[], values:double[])=
        if values = null then
            // set the structure of the jacobian
           // this particular jacobian is dense
            iRow.[0] <- 0
            jCol.[0] <- 0
            iRow.[1] <- 0
            jCol.[1] <- 1
            iRow.[2] <- 0
            jCol.[2] <- 2
            iRow.[3] <- 0
            jCol.[3] <- 3
            iRow.[4] <- 1
            jCol.[4] <- 0
            iRow.[5] <- 1
           jCol.[5] <- 1
           iRow.[6] <- 1
           jCol.[6] <- 2
           iRow.[7] <- 1
           jCol.[7] <- 3
       else
          // return the values of the jacobian of the constraints
          values.[0] <- x.[1] * x.[2] * x.[3] // 0,0
          values.[1] <- x.[0] * x.[2] * x.[3] // 0,1
          values.[2] <- x.[0] * x.[1] * x.[3] // 0,2
          values.[3] <- x.[0] * x.[1] * x.[2] // 0,3
          values.[4] <- 2.0 * x.[0] // 1,0
          values.[5] <- 2.0 * x.[1] // 1,1
          values.[6] <- 2.0 * x.[2] // 1,2
          values.[7] <- 2.0 * x.[3] // 1,3
       true

    member public this.eval_h (n:int, x:double[], new_x:bool, obj_factor:double,m:int, lambda:double[],  new_lambda:bool,nele_hess:int, iRow:int[], jCol:int[],values:double[]) =

        if values = null then
            // set the Hessian structure. This is a symmetric matrix, fill the lower left
            // triangle only.
            // the hessian for this problem is actually dense
            let mutable idx= 0; // nonzero element counter
            for row in 0..3 do
                 for col in 0..row do
                      iRow.[idx] <- row
                      jCol.[idx] <- col
                      idx<-idx+1
        else
            // return the values. This is a symmetric matrix, fill the lower left
            // triangle only
            // fill the objective portion
            values.[0] <- obj_factor * (2.0 * x.[3]); // 0,0
            values.[1] <- obj_factor * (x.[3]); // 1,0
            values.[2] <- 0.0; // 1,1
            values.[3] <- obj_factor * (x.[3]); // 2,0
            values.[4] <- 0.0; // 2,1
            values.[5] <- 0.0; // 2,2
            values.[6] <- obj_factor * (2.0 * x.[0] + x.[1] + x.[2]); // 3,0
            values.[7] <- obj_factor * (x.[0]); // 3,1
            values.[8] <- obj_factor * (x.[0]); // 3,2
            values.[9] <- 0.0; // 3,3
            // add the portion for the first constraint
            values.[1] <- values.[1]+ lambda.[0] * (x.[2] * x.[3]); // 1,0 */
            values.[3] <- values.[3]+ lambda.[0] * (x.[1] * x.[3]); // 2,0 */
            values.[4] <- values.[4]+ lambda.[0] * (x.[0] * x.[3]); // 2,1 */
            values.[6] <- values.[6]+ lambda.[0] * (x.[1] * x.[2]); // 3,0 */
            values.[7] <- values.[7]+ lambda.[0] * (x.[0] * x.[2]); // 3,1 */
            values.[8] <- values.[8]+ lambda.[0] * (x.[0] * x.[1]); // 3,2 */
            // add the portion for the second constraint
            values.[0] <- values.[0]+ lambda.[1] * 2.0; // 0,0 */
            values.[2] <- values.[2]+ lambda.[1] * 2.0; // 1,1 */
            values.[5] <- values.[5]+ lambda.[1] * 2.0; // 2,2 */
            values.[9] <- values.[9]+ lambda.[1] * 2.0; // 3,3 */
        true

// create the IpoptProblem

let p = new HS071()
// allocate space for the initial point and set the values
let x = [1.0; 5.0; 5.0; 1.0 ]
let mutable obj=0.0
let d_eval_g= new EvaluateConstraintsDelegate(fun (n:int) (x:double[]) (new_x:bool) (m:int) (g:double[])-> p.eval_g(n,x,new_x,m,g))
let d_eval_f= new EvaluateObjectiveDelegate(fun n x new_x obj_value-> p.eval_f(n,x,new_x,&obj_value))
let d_eval_grad_f= new EvaluateObjectiveGradientDelegate(fun n x new_x grad_f-> p.eval_grad_f(n,x,new_x,grad_f))
let d_eval_jac_g= new EvaluateJacobianDelegate(fun n x new_x m nele_jac iRow jCol values-> p.eval_jac_g(n,x,new_x,m,nele_jac,iRow,jCol,values))
let d_eval_h= new EvaluateHessianDelegate(fun n x new_x obj_factor m lambda new_Lambda nele_hess iRow jCol values-> p.eval_h(n,x,new_x,obj_factor,m,lambda,new_Lambda,nele_hess,iRow,jCol,values))

let problem = new Ipopt(p.n, p.x_L, p.x_U, p.m, p.g_L, p.g_U, p.nele_jac, p.nele_hess, d_eval_f, d_eval_g, d_eval_grad_f, d_eval_jac_g, d_eval_h)

// Set some options. Note the following ones are only examples, they might not be suitable for your problem.
problem.AddOption("tol", 1e-7) >ignore
problem.AddOption("mu_strategy", "adaptive") >ignore
problem.AddOption("output_file", "hs071_fs.txt") >ignore

// solve the problem
let status=problem.SolveProblem(x, &obj, null, null, null, null);
printfn "Optimization return status:%A" status
    for i in 0..3 do
        printfn "x[%d]=%f" i x.[i]
(*
Optimization return status:Solve_Succeeded
x[0]=1.000000
x[1]=4.743000
x[2]=3.821150
x[3]=1.379408
*)

Tuesday, September 14, 2010

Graphics in F# interactive

The following script snippets present graphics functionalities available in FlyingFrog.Graphics namespace from FSharpForVisualization assembly.
This functionalities allow the usage of F# interactive with plots much like other tools as Matlab or Python with MatPlotLib.
#I @"C:\cs\Docs"
#r "FSharp.PowerPack.dll"

// Reference WPF
#r "PresentationCore.dll"
#r "PresentationFramework.dll"
#r "WindowsBase.dll"

// Reference to FlyingFrog
#r @"FSharpForVisualization.dll"


open System.Windows
open System.Windows.Media
open FlyingFrog.Graphics;;

//Plot
Plot([Function sin], (-6., 6.));;

// Plot coordinates
let xys =
   [ for x in -5 .. 5 ->
       let x = float x
       x, sin x ]

Plot([Function sin; Data xys], (-6., 6.), Labels=(Text "x", Text "sin(x)"));;


// Progressive creation of a plot. Begin with an empty plot:
let g = Plot([], (-6., 6.));;

// And then specify a function to be plotted:
g.Data <- [Function sinc; Data[System.Math.PI, 0.]];;

// Specify as many functions:
g.Data <- [Function sin; Function cos];;

// Alter the styles:
g.Styles <- [Brushes.Black; Brushes.Gray];;



// 3D surface plot
let f x z =
    let r = 5. * sqrt(x*x + z*z)
    sin(r + 3. * x) / r

Plot3D(f, (-5., 5.), (-5., 5.), (-0.1, 1.));;
 
 
// Vector graphic shapes are styled by a sequence of strokes or fills:

let styles =

    [ Stroke(Brushes.Red, {width = 0.001}) ]
let geometry =

    let n = 50
    [ for i in 0 .. n do
         for j in 0 .. 1 do
            let x, y = float i / float n, float j
            yield Contour.Open(Point(x, y), [ Segment.Line(Point(y, 1. - x)) ]) ]
let scene = Shape(styles, geometry)

View scene;;



Plot([], (0., 1.), (0., 1.), Epilog=scene);;
 


// 3D object - icosahedron
let scene1 = Shape3D(Polyhedra.Icosahedron.mesh, Brushes.White)
View3D scene1;;

Sunday, September 12, 2010

Vector, Matrix and Linear Algebra in F#

The following script snippets present the vector and matrix functionalities in Microsoft.Fsharp.Numerics namespace from Fsharp.PowerPack assemble.
Also, they present the Linear Algebra functionalities from FlyingFrog namespace from FSharpForNumeric assemble.

//F# documentation at
//http://msdn.microsoft.com/library/dd233154(VS.100).aspx
//Fsharp.PowerPack documentation at
//http://research.microsoft.com/en-us/um/cambridge/projects/fsharp/manual/namespaces.html
//Flying frogs numerics documentation at
//http://www.ffconsultancy.com/products/fsharp_for_numerics/docs/namespaces.html


#I @"C:\Program Files (x86)\FSharpPowerPack-2.0.0.0\bin"
#r "Fsharp.PowerPack.dll"
#r "FSharpForNumerics.dll"
#r "FSharpForVisualization.dll"

open Microsoft.FSharp.Math
open FlyingFrog
open FlyingFrog.Graphics
open FlyingFrog.Graphics.Typeset;;

//vector (float)
let v1 = vector [1.0;1.0;1.0]
let v2 = vector [2.0;2.0;2.0]
let v3 = rowvec [3.0;3.0;3.0];;
val v1 : vector = vector [|1.0; 1.0; 1.0|]
val v2 : vector = vector [|2.0; 2.0; 2.0|]
val v3 : rowvec = rowvec [|3.0; 3.0; 3.0|]

//Simple operations
2.0*v1;;
v1+v2;;
v1-v2;;
v1.*v2;; //per element multiplication
let m1 = v2*v3;;

val it : Vector = vector [|2.0; 2.0; 2.0|]
val it : Vector = vector [|3.0; 3.0; 3.0|]
val it : Vector = vector [|-1.0; -1.0; -1.0|]
val it : Vector = vector [|2.0; 2.0; 2.0|]
val m1 : Matrix = matrix [[6.0; 6.0; 6.0]
[6.0; 6.0; 6.0]
[6.0; 6.0; 6.0]]
//Transpose
v1.Transpose;;
val it : RowVector = rowvec [|1.0; 1.0; 1.0|]

//Vector Generic
type intVec = Vector
let intVec ll = Vector.Generic.ofSeq ll
let A = intVec [1;2;3;4]
let B = intVec [1;1;1;1]
A+B;;
A*10;;
Vector.Generic.sum A;;

type intVec = Vector
val intVec : seq<'a> -> Vector<'a>
val A : Vector = vector [|1; 2; 3; 4|]
val B : Vector = vector [|1; 1; 1; 1|]

>val it : Vector = vector [|2; 3; 4; 5|]
> val it : Vector = vector [|10; 20; 30; 40|]
> val it : int = 10

//matrix (float)
let A = matrix [ [ 1.0; 7.0; 2.0 ];
[ 1.0; 3.0; 1.0 ];
[ 2.0; 9.0; 1.0 ]; ]
let B = matrix [ [ 10.0; 70.0; 20.0 ];
[ 10.0; 30.0; 10.0 ];
[ 20.0; 90.0; 10.0 ]; ];;

val A : matrix = matrix [[1.0; 7.0; 2.0]
[1.0; 3.0; 1.0]
[2.0; 9.0; 1.0]]
val B : matrix = matrix [[10.0; 70.0; 20.0]
[10.0; 30.0; 10.0]
[20.0; 90.0; 10.0]]
//print
let f (i:Matrix<_>)= printfn "%A" (i.ToArray2D())
let g (i:Vector<_>)= printfn "%A" i;;

val f : Matrix<'a> -> unit
val g : Vector<'a> -> unit

// sum, subtract
f (A+B)
> [[11.0; 77.0; 22.0]
[11.0; 33.0; 11.0]
[22.0; 99.0; 11.0]]
> [[-9.0; -63.0; -18.0]
[-9.0; -27.0; -9.0]
[-18.0; -81.0; -9.0]]

// matrix product
f (A*B);;

[[120.0; 460.0; 110.0]
[60.0; 250.0; 60.0]
[130.0; 500.0; 140.0]]

// element-wise product
f (A.*B);;

[[10.0; 490.0; 40.0]
[10.0; 90.0; 10.0]
[40.0; 810.0; 10.0]]
// scalar product
f (A * 2.0);;

[[2.0; 14.0; 4.0]
[2.0; 6.0; 2.0]
[4.0; 18.0; 2.0]]
f (2.0 * A);;
[[2.0; 14.0; 4.0]
[2.0; 6.0; 2.0]
[4.0; 18.0; 2.0]]
// negation of a matrix
f (-A);;

[[-1.0; -7.0; -2.0]
[-1.0; -3.0; -1.0]
[-2.0; -9.0; -1.0]]

// matrix-vector product
let b = vector [5.;8.;9.];;
g (A*b);;

val b : vector = vector [|5.0; 8.0; 9.0|]
> vector [|79.0; 38.0; 91.0|]

// Dimension
let dim = A.Dimensions;;

val dim : int * int = (3, 3)

//Transpose
Let A' = A.Transpose;;

val A' : Matrix = matrix [[1.0; 1.0; 2.0]
[7.0; 3.0; 9.0]
[2.0; 1.0; 1.0]]
//Number of rows,columns
let nrow = A.NumRows
let ncol = A.NumColslet;;

val nrow : int = 3
val ncol : int = 3

//Copy
let Anew = A.Copy();;

val Anew : Matrix = matrix [[1.0; 7.0; 2.0]
[1.0; 3.0; 1.0]
[2.0; 9.0; 1.0]]
// convert to a Array2D type
let Aarr = A.ToArray2D();;

val Aarr : float [,] = [[1.0; 7.0; 2.0]
[1.0; 3.0; 1.0]
[2.0; 9.0; 1.0]]
//Convert to vector (when the matrix has one row or column)
let Avec = A.ToVector()
let Arvec = A.ToRowVector();;

val Avec : Vector = vector [|1.0; 1.0; 2.0|]
val Arvec : RowVector = rowvec [|1.0; 7.0; 2.0|]

//operator [,]
A.[2,2];;
A.[2,2] <- 100.0;; A.[2,2];;

> val it : float = 100.0
> val it : unit = ()
> val it : float = 100.0

//sub-matrix
A.[1..2,1..2];;
A.[1..2,1..2] <- matrix [[2.;3.]; [8.;9.;]];; A.[1..2,1..2];;

> val it : Matrix = matrix [[3.0; 1.0]
[9.0; 1.0]]
> val it : unit = ()
> val it : Matrix = matrix [[2.0; 3.0]
[8.0; 9.0]]
//columns, rows
A.Column 2;;
A.Row 2;;
A.Columns (1,2);;
A.Rows (1,2);;

> val it : Vector = vector [|2.0; 3.0; 9.0|]
> val it : RowVector = rowvec [|2.0; 8.0; 9.0|]
> val it : Matrix = matrix [[7.0; 2.0]
[2.0; 3.0]
[8.0; 9.0]]
> val it : Matrix = matrix [[1.0; 2.0; 3.0]
[2.0; 8.0; 9.0]]
// sum and prod of all elements
let Asum = Matrix.sum A;;
let Aprod = Matrix.prod A;;

> val Asum : float = 35.0
> val Aprod : float = 12096.0

// create a matrix with 1
let C = Matrix.create 5 5 1.0;;

> val C : matrix = matrix [[1.0; 1.0; 1.0; 1.0; 1.0]
[1.0; 1.0; 1.0; 1.0; 1.0]
[1.0; 1.0; 1.0; 1.0; 1.0]
[1.0; 1.0; 1.0; 1.0; 1.0]
[1.0; 1.0; 1.0; 1.0; 1.0]]

// create a matrix with a function
let table = Matrix.init 3 3 (fun i j ->; (float i + 1.) * (float j + 1.));;

> val table : matrix = matrix [[1.0; 2.0; 3.0]
[2.0; 4.0; 6.0]
[3.0; 6.0; 9.0]]
// Identity
let I = Matrix.identity 2;;
val I : matrix = matrix [[1.0; 0.0]
[0.0; 1.0]]
// trace
let Atrace = Matrix.trace A;;

val Atrace : float = 12.0

//Map a function to the elements
let Asqr = Matrix.map (fun x -> x*x) A;;

> val Asqr : matrix = matrix [[1.0; 49.0; 4.0]
[1.0; 4.0; 9.0]
[4.0; 64.0; 81.0]]
//Sparse matrix
let D = Matrix.initSparse 3 3 [ (0,0,1.0); (1,1,2.0); (2,2,3.0); ];;

> val D : matrix = matrix [[1.0; 0.0; 0.0]
[0.0; 2.0; 0.0]
[0.0; 0.0; 3.0]]
D.IsSparse;;
> val it : bool = true

//Matrix Generic
type intmatrix = Matrix
let intmatrix ll = Matrix.Generic.ofSeq ll
let C = intmatrix [ [1;2]; [3;4] ]
let D = intmatrix [ [1;1]; [1;1] ]
Matrix.Generic.inplaceAdd C D;;
C+D;;
C * 10;;
type intmatrix = Matrix
val intmatrix : seq<#seq<'b>> -> Matrix<'b>
val C : Matrix = matrix [[2; 3]
[4; 5]]
val D : Matrix = matrix [[1; 1]
[1; 1]]
> val it : Matrix = matrix [[3; 4]
[5; 6]]
> val it : Matrix = matrix [[20; 30]
[40; 50]]
//Linear Algebra - float

//Solve linear equations
let A1 =
[ [1; 2; 3];
[1; 4; 9];
[1; 8; 27] ]
let A2 = Matrix.ofSeq(Seq.map (Seq.map float) A1)
let b = vector [2.; 3.; 4.]
let x = LinearAlgebra.Float.linear_solve A b;;

val A1 : int list list = [[1; 2; 3]; [1; 4; 9]; [1; 8; 27]]
val A2 : matrix = matrix [[1.0; 2.0; 3.0]
[1.0; 4.0; 9.0]
[1.0; 8.0; 27.0]]
val b : vector = vector [|2.0; 3.0; 4.0|]
val x : vector = vector [|0.5; 1.0; -0.1666666667|]

//Determinant
LinearAlgebra.Float.determinant A2;;

> val it : float = -12.0

//Inverse
LinearAlgebra.Float.inverse A2;;
> val it : matrix = matrix [[3.0; -2.5; 0.5]
[-1.5; 2.0; -0.5]
[0.3333333333; -0.5; 0.1666666667]]
//svd
LinearAlgebra.Float.svd A2;;

> val it : matrix * matrix * matrix =
(matrix [[0.1166955507; 0.6917676817; 0.7126286712]
[0.3269819117; 0.6507677484; -0.6852621157]
[0.9377979409; -0.3129837252; 0.1502538182]],
matrix [[30.0404326; 0.0; 0.0]
[0.0; 1.878075836; 0.0]
[0.0; 0.0; 0.2126972812]],
matrix [[0.0459872007; 0.5481949583; 0.835085304]
[0.301051; 0.7894977011; -0.5348473383]
[0.9524985421; -0.2759993977; 0.1287278512]])

//Linear Algebra – rational

//Solve linear equations
let A3 =
[ [1; 2; 3];
[1; 4; 9];
[1; 8; 27] ]
let A4 = Matrix.Generic.ofSeq(Seq.map (Seq.map BigNum.FromInt) A3)
let b1 = Vector.Generic.ofSeq [2N; 3N; 4N]
let x1 = LinearAlgebra.Rational.linear_solve A4 b1;;
A4 * x1;;

> val A3 : int list list = [[1; 2; 3]; [1; 4; 9]; [1; 8; 27]]
val A4 : Matrix = matrix [[1N; 2N; 3N]
[1N; 4N; 9N]
[1N; 8N; 27N]]
val b1 : Vector = vector [|2N; 3N; 4N|]
val x1 : Vector = vector [|1/2N; 1N; -1/6N|]
> val it : Vector = vector [|2N; 3N; 4N|]


//Inverse - rational
LinearAlgebra.Rational.inverse A4;;

> val it : Matrix = matrix [[3N; -5/2N; 1/2N]
[-3/2N; 2N; -1/2N]
[1/3N; -1/2N; 1/6N]]

//Linear Algebra - complex

type cmatrix = Matrix
let cmatrix ll = Matrix.Generic.ofSeq ll

let A5 = cmatrix [ [ complex 1.0 1.0; complex 1.0 1.0];
[ complex 1.0 1.0; complex 1.0 1.0]];;

> type cmatrix = Matrix
val cmatrix : seq<#seq<'b>> -> Matrix<'b>
val A5 : Matrix = matrix [[1r+1i; 1r+1i]
[1r+1i; 1r+1i]]
//Determinant - complex
LinearAlgebra.Complex.determinant A5;;

> val it : Complex = 0r+0i {Conjugate = 0r+0i;
ImaginaryPart = 0.0;
Magnitude = 0.0;
Phase = 0.0;
RealPart = 0.0;
i = 0.0;
r = 0.0;}
f (A-B);;

Saturday, September 11, 2010

F#- Introduction

F# is a first class language in .NET and Visual Studio together with C#, VB and C++. Solutions are created and debugged in the same way. IntelliSense is available. Generated assembles can be used in the .NET solutions. Execution speed is equivalent to C#.

But, F# also has unique features as F# Interactive and F# Scripting Files. These features are attractive to scientific and engineering developers that are used to tools like Matlab or Python.

Installation

A free installation is available as follow:

1) Install .NET 4.0 from
http://www.microsoft.com/downloads/en/details.aspx?FamilyID=9cfb2d51-5ff4-4491-b0e5-b386f32c0992&displaylang=en

2) Install Visual Studio Shell 10 from

http://www.microsoft.com/downloads/en/details.aspx?familyid=8E5AA7B6-8436-43F0-B778-00C3BCA733D3&displaylang=en

3) Install F# 2.0 from

http://www.microsoft.com/downloads/en/details.aspx?FamilyID=f8c623ae-aef6-4a06-a185-05f59be47d67&displaylang=en

4) Add extra content from:

http://fsharppowerpack.codeplex.com/releases/view/40168

Using F# Interactive

In order to open an interactive session, start the following executable:

C:\Program Files (x86)\Microsoft F#\v4.0\fsi.exe

Interactive session is also available inside Visual Studio. You can be shown it with ‘View/Other Windows/F# Interactive’.

Using F# Script Files

Script files have extension ‘.fsx’. You can create them with Notepad. By right-clicking on Windows Explorer, you can execute it by selecting “Run with F# Interactive”. This will start the script file with FSI.

Using Visual Studio to create script files has the advantage of IntelliSense. Create an script by using ‘File/New/File/F# Script File’.

Compiling F#

In order to generate an executable or a library, create a ‘.fs’ file in Notepad and use it with:

C:\Program Files\Microsoft F#\v4.0\fsc.exe

Using Visual Studio to create source files has the advantage of IntelliSense and easy reference to other assembles and .NET resources. Create a project by using ‘File/New/Project/F# Application or F# Library’.

Using Source File and Interactive

Code can be entered in Visual Studio source file. You can transfer and execute sections of the code by marking it and pressing ‘Alt+Enter’. This is a method to debug the code, traditionally used by systems with interactive capability.