Change over to gsl-ocaml

lib/CamlPaml/Fit.ml

@@ -10,13 +10,13 @@ class maximizer f init lo hi =
 			| any ->
 				exn := Some any
 				nan
-	let minimizer = Gsl_min.make Gsl_min.BRENT f' ~min:init ~lo:lo ~up:hi
+	let minimizer = Gsl.Min.make Gsl.Min.BRENT f' ~min:init ~lo:lo ~up:hi
 
 	object
-		method maximum () = Gsl_min.minimum minimizer
-		method interval () = Gsl_min.interval minimizer
+		method maximum () = Gsl.Min.minimum minimizer
+		method interval () = Gsl.Min.interval minimizer
 		method iterate () =
-			Gsl_min.iterate minimizer
+			Gsl.Min.iterate minimizer
 			match !exn with
 				| Some ex -> raise ex
 				| None -> ()
@@ -55,7 +55,7 @@ class multi_maximizer f df init =
 	let f' ~x =
 		if !exn = None then
 			try
-				0. -. f (Gsl_vector.to_array x)
+				0. -. f (Gsl.Vector.to_array x)
 			with
 				| any ->
 					exn := Some any
@@ -65,31 +65,31 @@ class multi_maximizer f df init =
 	let df' ~x ~g =
 		if !exn = None then
 			try
-				let dfv = Gsl_vector.of_array (df (Gsl_vector.to_array x))
-				Gsl_vector.scale dfv (-1.) 
-				Gsl_vector.set_zero g
-				Gsl_vector.add g dfv
+				let dfv = Gsl.Vector.of_array (df (Gsl.Vector.to_array x))
+				Gsl.Vector.scale dfv (-1.) 
+				Gsl.Vector.set_zero g
+				Gsl.Vector.add g dfv
 			with
 				| any -> exn := Some any
 	let gsl_fdf = {
-		Gsl_fun.multim_f = f';
-		Gsl_fun.multim_df = df';
-		Gsl_fun.multim_fdf = (fun ~x ~g -> let rslt = f' ~x in (if !exn = None then (df' ~x ~g; rslt) else rslt));
+		Gsl.Fun.multim_f = f';
+		Gsl.Fun.multim_df = df';
+		Gsl.Fun.multim_fdf = (fun ~x ~g -> let rslt = f' ~x in (if !exn = None then (df' ~x ~g; rslt) else rslt));
 	}
 	
 	
-	let minimizer = Gsl_multimin.Deriv.make Gsl_multimin.Deriv.VECTOR_BFGS2 n gsl_fdf ~x:(Gsl_vector.of_array init) ~step:1. ~tol:1e-6
+	let minimizer = Gsl.Multimin.Deriv.make Gsl.Multimin.Deriv.VECTOR_BFGS2 n gsl_fdf ~x:(Gsl.Vector.of_array init) ~step:1. ~tol:1e-6
 	
 	object
 		method maximum () = 
-			let x = Gsl_vector.create n
-			let g = Gsl_vector.create n
-			let opt = 0. -. Gsl_multimin.Deriv.minimum ~x:x ~g:g minimizer
-			Gsl_vector.scale g (-1.)
+			let x = Gsl.Vector.create n
+			let g = Gsl.Vector.create n
+			let opt = 0. -. Gsl.Multimin.Deriv.minimum ~x:x ~g:g minimizer
+			Gsl.Vector.scale g (-1.)
 
-			opt, (Gsl_vector.to_array x), (Gsl_vector.to_array g)
+			opt, (Gsl.Vector.to_array x), (Gsl.Vector.to_array g)
 		method iterate () =
-			Gsl_multimin.Deriv.iterate minimizer
+			Gsl.Multimin.Deriv.iterate minimizer
 			match !exn with Some ex -> raise ex | None -> ()
 			
 let make_multi_maximizer ~f ~df ~init = (new multi_maximizer f df init)

lib/CamlPaml/P.ml

@@ -1,15 +1,15 @@
 open Printf
 
-type matrix = Gsl_matrix.matrix
+type matrix = Gsl.Matrix.matrix
 
 let validate ?(tol=1e-6) p =
-	let (n,n') = Gsl_matrix.dims p
+	let (n,n') = Gsl.Matrix.dims p
 	
 	if n <> n' then invalid_arg "CamlPaml.P.validate: non-square matrix"
 	if n < 2 then invalid_arg "CamlPaml.P.validate: trivial matrix"
 	
 	for i = 0 to n-1 do
-		let pi = Gsl_matrix.row p i
+		let pi = Gsl.Matrix.row p i
 		let rowsum = ref 0.
 		for j = 0 to n-1 do
 			if pi.{j} < 0. || pi.{j} > 1. then invalid_arg (sprintf "CamlPaml.P.validate: P[%d,%d] = %f" i j pi.{j})

lib/CamlPaml/P.mli

@@ -1,6 +1,6 @@
 (** substitution probability matrices (P matrices) *)
 
-type matrix =  (float, Bigarray.float64_elt, Bigarray.c_layout) Bigarray.Array2.t (** [Gsl_matrix.matrix] *)
+type matrix =  (float, Bigarray.float64_elt, Bigarray.c_layout) Bigarray.Array2.t (** [Gsl.Matrix.matrix] *)
 
 (** validate that given array is a positive square matrix in which the rows sum to 1
 @raise Invalid_argument if not*)

lib/CamlPaml/PhyloEM.ml

@@ -37,7 +37,7 @@ let branch_ell m ss br =
 	let pbr = PM.p m br
 
 	for k = 0 to n-1 do
-		let pbrk = Gsl_matrix.row pbr k
+		let pbrk = Gsl.Matrix.row pbr k
 		let ssbrk = ssbr.(k)
 		for l = 0 to n-1 do
 			let ssbrkl = ssbrk.(l)
@@ -127,7 +127,7 @@ let d_ell_dQ_dxi inst ss i =
 			for a = 0 to n-1 do
 				let ssbra = ssbr.(a)
 				let dPt_dxi_a = dPt_dxi.(a)
-				let pa = Gsl_matrix.row p a
+				let pa = Gsl.Matrix.row p a
 				for b = 0 to n-1 do
 					let ssbrab = ssbra.(b)
 					if ssbrab > 0. then
@@ -156,7 +156,7 @@ let d_ell_dbranch inst br ss =
 		for a = 0 to n-1 do
 			let ssbra = ssbr.(a)
 			let dPt_dt_a = dPt_dt.(a)
-			let pa = Gsl_matrix.row p a
+			let pa = Gsl.Matrix.row p a
 			for b = 0 to n-1 do
 				let ssbrab = ssbra.(b)
 				if ssbrab > 0. then

lib/CamlPaml/PhyloLik.ml

@@ -8,39 +8,39 @@ module IntSet = Set.Make(Int)
 
 type leaf = [`Certain of int | `Distribution of float array | `Marginalize]
 
-let raw_leaf (k,i) = Gsl_vector.of_array (Array.init k (fun j -> if i = j then 1. else 0.))
-let raw_marg k = Gsl_vector.of_array (Array.make k 1.)
+let raw_leaf (k,i) = Gsl.Vector.of_array (Array.init k (fun j -> if i = j then 1. else 0.))
+let raw_marg k = Gsl.Vector.of_array (Array.make k 1.)
 let hashcons_raw_leaf = Tools.weakly_memoize raw_leaf
 let hashcons_raw_marg = Tools.weakly_memoize raw_marg
 
 let get_leaf k leaf = match leaf with
 	| `Certain i -> hashcons_raw_leaf (k,i)
-	| `Distribution ar -> Gsl_vector.of_array ar
+	| `Distribution ar -> Gsl.Vector.of_array ar
 	| `Marginalize -> hashcons_raw_marg k
 
 type workspace = {
 	mutable generation : int;
-	data : Gsl_matrix.matrix
+	data : Gsl.Matrix.matrix
 }
 
 type intermediate = {
 	tree : T.t;
-	pms : Gsl_matrix.matrix array;
+	pms : Gsl.Matrix.matrix array;
 	leaves : leaf array;
 	
 	workspace : workspace;
 	my_generation : int;
 	
-	alpha : Gsl_matrix.matrix; (** actually a sub-matrix of workspace.data *)
+	alpha : Gsl.Matrix.matrix; (** actually a sub-matrix of workspace.data *)
 	mutable z : float;
 	mutable have_alpha : bool;
-	beta : Gsl_matrix.matrix;  (** actually a sub-matrix of workspace.data *)
+	beta : Gsl.Matrix.matrix;  (** actually a sub-matrix of workspace.data *)
 	mutable have_beta : bool;
 }
 
 let new_workspace tree dim =
 	let rows = 2 * (T.size tree) - (T.leaves tree)
-	{ generation = min_int; data = Gsl_matrix.create rows dim }
+	{ generation = min_int; data = Gsl.Matrix.create rows dim }
 
 let empty = [||]
 let prepare ?workspace tree pms prior leaves =
@@ -54,7 +54,7 @@ let prepare ?workspace tree pms prior leaves =
 	let workspace = match workspace with Some x -> x | None -> new_workspace tree k
 	workspace.generation <- (if workspace.generation = max_int then min_int else workspace.generation+1)
 	
-	let rows,cols = Gsl_matrix.dims workspace.data
+	let rows,cols = Gsl.Matrix.dims workspace.data
 	if rows < (2*n-nl) || cols <> k then invalid_arg "CamlPaml.Infer.prepare: inappropriate workspace dimensions"
 	let alpha = Bigarray.Array2.sub_left workspace.data 0 (n-nl)
 	let beta = Bigarray.Array2.sub_left workspace.data (n-nl) n
@@ -66,15 +66,15 @@ let prepare ?workspace tree pms prior leaves =
 
 (* Inside algo (aka Felsenstein algo.) *)
 let alpha_get x br =
-	let k = snd (Gsl_matrix.dims x.alpha)
+	let k = snd (Gsl.Matrix.dims x.alpha)
 	let nl = T.leaves x.tree
-	if br < nl then get_leaf k x.leaves.(br) else Gsl_matrix.row x.alpha (br-nl)
+	if br < nl then get_leaf k x.leaves.(br) else Gsl.Matrix.row x.alpha (br-nl)
 	
 let ensure_alpha x =
 	if not x.have_alpha then
 		let n = T.size x.tree
 		let nl = T.leaves x.tree
-		let k = snd (Gsl_matrix.dims x.alpha)
+		let k = snd (Gsl.Matrix.dims x.alpha)
 		for i = nl to n-1 do
 			let (lc,rc) = T.children x.tree i
 			assert (lc >= 0 && rc >= 0)
@@ -85,35 +85,35 @@ let ensure_alpha x =
 			let arc = alpha_get x rc
 
 			for a = 0 to k-1 do
-				let lsa = Gsl_matrix.row ls a
-				let rsa = Gsl_matrix.row rs a
-				x.alpha.{i-nl,a} <- (Gsl_blas.dot lsa alc) *. (Gsl_blas.dot rsa arc)
+				let lsa = Gsl.Matrix.row ls a
+				let rsa = Gsl.Matrix.row rs a
+				x.alpha.{i-nl,a} <- (Gsl.Blas.dot lsa alc) *. (Gsl.Blas.dot rsa arc)
 
-		x.z <- Gsl_blas.dot (alpha_get x (n-1)) (Gsl_matrix.row x.beta (n-1))
+		x.z <- Gsl.Blas.dot (alpha_get x (n-1)) (Gsl.Matrix.row x.beta (n-1))
 		x.have_alpha <- true
 
 (* Outside algo *)
 let ensure_beta x =
 	ensure_alpha x
 	if not x.have_beta then
-		let k = snd (Gsl_matrix.dims x.beta)
-		let inter = Gsl_vector.create k
-		let ps_colb = Gsl_vector.create k
+		let k = snd (Gsl.Matrix.dims x.beta)
+		let inter = Gsl.Vector.create k
+		let ps_colb = Gsl.Vector.create k
 		for i = (T.root x.tree)-1 downto 0 do
 			let p = T.parent x.tree i
 			assert (p > i)
 			let s = T.sibling x.tree i
 			let ps = x.pms.(i)
 			let ss = x.pms.(s)
-			let bp = Gsl_matrix.row x.beta p
+			let bp = Gsl.Matrix.row x.beta p
 			let xas = alpha_get x s
 
 			for a = 0 to k-1 do
-				inter.{a} <- bp.{a} *. (Gsl_blas.dot (Gsl_matrix.row ss a) xas)
+				inter.{a} <- bp.{a} *. (Gsl.Blas.dot (Gsl.Matrix.row ss a) xas)
 
 			for b = 0 to k-1 do
 				for a = 0 to k-1 do ps_colb.{a} <- ps.{a,b}
-				x.beta.{i,b} <- Gsl_blas.dot inter ps_colb
+				x.beta.{i,b} <- Gsl.Blas.dot inter ps_colb
 		x.have_beta <- true
 		
 let likelihood x =
@@ -124,21 +124,21 @@ let likelihood x =
 let node_posterior inferred i =
 	if inferred.workspace.generation <> inferred.my_generation then failwith "CamlPaml.PhyloLik.node_posterior: invalidated workspace"
 	ensure_beta inferred
-	let k = snd (Gsl_matrix.dims inferred.alpha)
+	let k = snd (Gsl.Matrix.dims inferred.alpha)
 	if inferred.z = 0. then
 		Array.make k 0.
 	else
 		let nleaves = T.leaves inferred.tree
 
 		if i < nleaves then
-			Gsl_vector.to_array (alpha_get inferred i)
+			Gsl.Vector.to_array (alpha_get inferred i)
 		else
 			Array.init k (fun x -> (alpha_get inferred i).{x} *. inferred.beta.{i,x} /. inferred.z)
 
 let add_branch_posteriors ?(weight=1.0) inferred branch ecounts =
 	if inferred.workspace.generation <> inferred.my_generation then failwith "CamlPaml.PhyloLik.add_branch_posterior: invalidated workspace"
 	ensure_beta inferred
-	let k = snd (Gsl_matrix.dims inferred.alpha)
+	let k = snd (Gsl.Matrix.dims inferred.alpha)
 	
 	if Array.length ecounts <> k then invalid_arg "CamlPaml.Infer.add_branch_posteriors: wrong matrix dimension"
 	if branch < 0 || branch >= (T.root inferred.tree) then invalid_arg "CamlPaml.Infer.add_branch_posteriors: invalid branch"
@@ -150,7 +150,7 @@ let add_branch_posteriors ?(weight=1.0) inferred branch ecounts =
 		let p = T.parent tree branch
 		let sib = T.sibling tree branch
 
-		let bp = Gsl_matrix.row inferred.beta p
+		let bp = Gsl.Matrix.row inferred.beta p
 		let ab = alpha_get inferred branch
 
 		let xas = alpha_get inferred sib
@@ -159,10 +159,10 @@ let add_branch_posteriors ?(weight=1.0) inferred branch ecounts =
 		for a = 0 to k-1 do
 			let bpa = bp.{a}
 			if bpa > 0. then
-				let sma = Gsl_matrix.row sm a
-				let smsa = Gsl_matrix.row sms a
+				let sma = Gsl.Matrix.row sm a
+				let smsa = Gsl.Matrix.row sms a
 
-				let bpa_sibtot = bpa *. (Gsl_blas.dot smsa xas)
+				let bpa_sibtot = bpa *. (Gsl.Blas.dot smsa xas)
 
 				let ecountsa = ecounts.(a)
 				if Array.length ecountsa <> k then invalid_arg "CamlPaml.Infer.add_branch_posteriors: wrong matrix dimension"
@@ -172,7 +172,7 @@ let add_branch_posteriors ?(weight=1.0) inferred branch ecounts =
 					ecountsa.(b) <- ecountsa.(b) +. weight *. pr
 
 let branch_posteriors inferred branch =
-	let k = snd (Gsl_matrix.dims inferred.alpha)
+	let k = snd (Gsl.Matrix.dims inferred.alpha)
 	let a = Array.make_matrix k k 0.
 	add_branch_posteriors inferred branch a
 	a

lib/CamlPaml/PhyloModel.ml

@@ -33,7 +33,7 @@ let prepare_lik ?workspace m leaves = PhyloLik.prepare ?workspace m.tree m.pms m
 let checksum = 1., 1e-6
 
 let simulate m =
-	let branch_choosers = m.pms |> Array.map (fun pm -> (Array.map (Tools.random_chooser ~checksum:checksum) (Gsl_matrix.to_arrays pm)))
+	let branch_choosers = m.pms |> Array.map (fun pm -> (Array.map (Tools.random_chooser ~checksum:checksum) (Gsl.Matrix.to_arrays pm)))
 	let root_chooser = Tools.random_chooser ~checksum:checksum m.prior
 	fun ?root ?a () ->
 		let t = m.tree