hubFS: THE place for F#

Hello,
I've worked some time ago in a bridge of F# and the dnAnalytics numerical library (see previous posts on this subject).

Yesterday I revamped the code so it is now 0.2 compliant (http://www.dnanalytics.net/files/2/default.aspx)
You'll find attached , and below, the library .
The goal is to provide routines and scripting in the 'matlab-like way of life'.
You are also able to add a matrix and a scalar, to do element by element operations, etc...

As an example, I took a full-fledged ordinary least square esimtation, so You'll have 2 examples of OLS with Marcus example in http://www.dnanalytics.net/blogs/marcus/archive/2006/05/07/189.aspx

Feel free to comment, and criticize.

(Sorry I haven't been able to attach any file for the moment)


Regards

Julien Laugel



An example :
#r "dnAnalytics.dll";;

#r "matrixlib.dll";;


open Matrixlib

open dnAnalytics.LinearAlgebra.IO


//Take data

let a = new DelimitedMatrixReader(";");;

let m = M(a.ReadMatrix("c:\mat1.csv"));;

//1st col = input

let y = submat m 0 156 0 0 ;;

//last cols = regressors

let x = submat m 0 156 1 3 ;;

//Result : ols_res

//See library to view all kind of output

let ols_res = ols y x ;;

ols_res.Beta;;

ols_res.Rbar;;

ols_res.Resid;;


Here is the library :
module Matrixlib

//#r "dnAnalytics.dll";;


open dnAnalytics.LinearAlgebra;;

open dnAnalytics.LinearAlgebra.Decomposition;;



open System;;

(*This is the type everywhere in the library*)

type genmat = M of Matrix


(*genmat -> Array conversion*)

let g2a (M a ) = a.ToArray();;

(*genmat -> double conversion (for scalar) *)

let g2s (M a ) =

let b = g2a (M a) in

b.(0,0);;


(* init a matrix with a function of its coordinates *)

(*Exampe : let m1 = init 5 5 (fun i j -> float(i) + float(j) );;*)

let init m n (f : int -> int -> float ) =

let tmp = Array2.init m n f in

let mat_tmp = MatrixBuilder.CreateMatrix(m,n) in

let _ = Array2.iteri (fun i j x -> (mat_tmp.Item(i,j)<- x) ) tmp in

M(mat_tmp) ;;


(*Declare a scalar (e.g. 1.) as a genmat*)

let scamat (a : float ) =

let m = MatrixBuilder.CreateMatrix(1,1) in

let _ = m.Item(0,0) <- a in

M(m);;


(*Declares a float array as a genmat

It is a (n,1) matrix *)

let vecmat (a : float array) =

let m = MatrixBuilder.CreateMatrix(a.Length,1) in

for i=0 to (a.Length - 1) do

m.Item(i,0) <- a.(i)

done;

M(m);;


(*Vector to Matrix conversion

* In some place in dNA, you have to use vectors

* In matlab-like languages, a vector is a special case of a matrix

* *)


let vector_to_mat (a:Vector) =

let m = MatrixBuilder.CreateMatrix(a.Count,1) in

for i=0 to (a.Count - 1) do

m.Item(i,0) <- a.Item(i)

done;

M(m);;


(*Boolean function :

* ckeck the shape of the genmat*)

let is_vec (a:Matrix) = (a.Columns = 1);; //column vector?

let is_rowvec (a:Matrix) = (a.Rows = 1);; //row vector?

let is_sca (a:Matrix) = ((a.Rows = 1) & (a.Columns = 1));; //scalar ?


(*vector indexing : you just need one dimension

* Uses the new method SubMatrix

TODO : boolean indexing

uses Range objects : init->end *)

let subvec (M a ) i j =

if (is_vec a) then M(a.SubMatrix(new Range(i,j),new Range(0,0)))

else if is_rowvec a then M(a.SubMatrix(new Range(0,0),new Range(i,j)))

else failwith "not implemented yet" ;;

(*An easier way to do vector indexing :

* e.g. my_vector.[0,5]*)

let ( .[] ) m (i,j) = subvec m i j;;


(*Submatrix extraction*)

(* indexing order is : beginning row, ending row, beginning column, ending column *)

let submat (M a) i1 i2 j1 j2 = M(a.SubMatrix(new Range(i1,i2),new Range(j1,j2)));;



(*Basic Operations

* You can use it with different couples of matrix and scalar *)


let add (M a) (M b) =

if ((a.Columns = b.Columns) & (a.Rows = b.Rows)) then

M(Matrix.op_Addition(a,b))

else if is_sca b then

let m = MatrixBuilder.CreateMatrix(a.Rows,a.Columns) in

let _ = a.Add(b.Item(0,0),m) in

M(m)

else if is_sca a then

let m = MatrixBuilder.CreateMatrix(b.Rows,b.Columns) in

let _ = b.Add(a.Item(0,0),m) in

M(m)

else failwith "error"


let sub (M a) (M b) =

if ((a.Columns = b.Columns) & (a.Rows = b.Rows)) then

M(Matrix.op_Subtraction(a,b))

else if is_sca b then

let m = MatrixBuilder.CreateMatrix(a.Rows,a.Columns) in

let _ = a.Add(-b.Item(0,0),m) in

M(m)

else if is_sca a then

let m = MatrixBuilder.CreateMatrix(b.Rows,b.Columns) in

let _ = b.Add(-a.Item(0,0),m) in

M(m)

else failwith "error"


let mult (M a) (M b) =

if is_sca b then M(Matrix.op_Multiply(a,b.Item(0,0)))

else if is_sca a then M(Matrix.op_Multiply(a.Item(0,0),b))

else M(Matrix.op_Multiply(a,b))


let divide (M a) (M b) =

if is_sca b then M(Matrix.op_Division(a,b.Item(0,0)))

else failwith "error"



(*Classical methods *)

let rows (M a) = a.Rows

let cols (M a) = a.Columns

let conformable a b = (rows a = rows b) & (cols a = cols b) //are the two matrices the same size ?



let zeros m n = M(MatrixBuilder.CreateMatrix( m , n , 0.));; //Build a matrix of zeros

let ones m n = M(MatrixBuilder.CreateMatrix( m , n , 1.));; //Build a matrix of ones


(*Mapping functions*)


(*Apply a float function to each matrix elt*)

let map (f : float -> float) (M x) =

let y = x.Clone() in

for i = 0 to (x.Rows -1) do

for j = 0 to (x.Columns -1) do

y.Item(i,j) <- f (x.Item(i,j))

done;

done;

M y;;


(*If 2 matrices are of same size, apply a function taking as args each elts of

* the 2 matrices

* kinda general elt to elt matrix operation *)

let map2 (f : float -> float -> float) (M a) (M b) =

if conformable (M a) (M b) then

let z = a.Clone() in

for i =0 to (a.Rows -1) do

for j = 0 to (a.Columns -1) do

z.Item(i , j) <- f (a.Item(i,j)) (b.Item(i,j))

done;

done;

M z;

else failwith "Matrices not conformable"


(*Another map function : the output is an 'a array*)

let map_to_array (f : float -> 'a) (M x) =

let y = Array.create (x.Rows) (f (x.Item(0,0))) in

for i = 0 to (x.Rows -1) do

for j = 0 to (x.Columns -1) do

y.(i) <- f (x.Item(i,j))

done;

done;

y;;



(*dnA library wrappings*)

let qr_decomp (M a) = new QR(a)

let qr_Q (a:QR) = M(a.Q())

let qr_R (a:QR) = M(a.R())

let qr_solve (a:QR) (M b) = M(a.Solve(b))


(*Common matrices*)


let eye n = let tmp = MatrixBuilder.CreateIdentityMatrix(n) in

M(tmp);;


let diag (M a) = vector_to_mat(a.Diagonal())

//rand is a uniform random

let rand m n =

let rng = new Random() in init m n (fun i j -> rng.NextDouble()) ;;


(*common matrix operations*)

let inverse (M a) = M(a.Inverse()) //inverting a matrix

let mldivide (a) (b) = mult (inverse a) b // the backslash matlab-esque operator

let (.^) a (x:float) = map (fun t -> t ** x) a //Power each elt matrix


(*Sum each column of the genmat*)

let sum (M a) =

let m = MatrixBuilder.CreateMatrix(1 , a.Columns , 0.) in

for j = 0 to (a.Columns -1) do

for i = 0 to (a.Rows -1) do

m.Item(0,j) <- m.Item(0,j) + a.Item(i,j)

done;

done;

M(m);;

(*Mean of each column of the genmat*)

let mean (a) =

let s = sum a in divide s (scamat (float (rows a)))

(* "Dotted operators" :

* They are classical operators operating on e element by element basis *)


let dot_add = map2 (+)

let dot_minus = map2 (-)

let dot_slash = map2 (/)

let dot_times = map2 ( * )

let transpose (M a) = M(a.Transpose())

let (~!) a = transpose a //Yet another way to transpose

type genmat

with

static member ( + ) ((m1:genmat),(m2:genmat)) = add m1 m2

static member ( - ) ((m1:genmat),(m2:genmat)) = sub m1 m2

static member ( * ) ((m1:genmat),(m2:genmat)) = mult m1 m2

static member ( / ) ((m1:genmat),(m2:genmat)) = divide m1 m2

static member ( .+ ) ((m1:genmat),(m2:genmat)) = dot_add m1 m2

static member ( .- ) ((m1:genmat),(m2:genmat)) = dot_minus m1 m2

static member ( .* ) ((m1:genmat),(m2:genmat)) = dot_times m1 m2

static member ( ./ ) ((m1:genmat),(m2:genmat)) = dot_slash m1 m2

end


(*Example : building a full-fledged ordinary least square function,

using genmat. You have residuals, correlations, betas, significance values,

etc...*)

type ols = {Beta : genmat ; Resid : genmat ; Sige : genmat ; Yhat :genmat ; Tstat : genmat ; Rbar : genmat ; Rsqr : genmat ; Dw :genmat};;


(*Ordinary Least Square regression *)

(*I took the code from Jim Lesage cool Econometrics ToolBox

* All potential errors are mine

* the code only uses genmat type*)


let ols (y) (x) =

let nb_obs = rows x in

let nb_var = cols x in

let reg = qr_decomp x in

let q = qr_Q reg in

let beta = qr_solve reg y in

let r = qr_R reg in

let yhat = x * beta in

let tmp_ = (~!r) * r in

let id = eye(nb_var) in

let xpxi = mldivide tmp_ id in

let resid = y - yhat in

let sigu = (~! resid) * resid in

let sige = (sigu * scamat(1./float_of_int(nb_obs - nb_var))) in

let tmp = sige * (diag xpxi) in

let tstat = beta ./ (tmp .^ 0.5) in

let ym = y - (mean y) in

let rsqr1 = sigu in

let rsqr2 = (~! ym) * ym in

let rsqr = (scamat 1.0) - (rsqr1 / rsqr2) in

let rsqr1 = rsqr1 / scamat ((float_of_int nb_obs) - (float_of_int nb_var)) in

let rsqr2 = rsqr2 / scamat ((float_of_int nb_obs) - 1.0) in

let rbar = (scamat 1.) - (rsqr1 / rsqr2) in

let ediff = (subvec resid 1 (nb_obs - 1)) - (subvec resid 0 (nb_obs - 2)) in

let dw = (~! ediff) * ediff in

{ Beta = beta ; Resid = resid ; Sige = sige ; Yhat = yhat ; Tstat =tstat ; Rbar = rbar ; Rsqr = rsqr ; Dw = dw };;

Hi nice work.

I have problems with this part:


type genmat
with
static member ( + ) ((m1:genmat),(m2:genmat)) = add m1 m2
static member ( - ) ((m1:genmat),(m2:genmat)) = sub m1 m2
static member ( * ) ((m1:genmat),(m2:genmat)) = mult m1 m2
static member ( / ) ((m1:genmat),(m2:genmat)) = divide m1 m2
static member ( .+ ) ((m1:genmat),(m2:genmat)) = dot_add m1 m2
static member ( .- ) ((m1:genmat),(m2:genmat)) = dot_minus m1 m2
static member ( .* ) ((m1:genmat),(m2:genmat)) = dot_times m1 m2
static member ( ./ ) ((m1:genmat),(m2:genmat)) = dot_slash m1 m2
end

Any suggestions?