From 9d8fa1f94b2d000981d44d616a0371c992f9a9d4 Mon Sep 17 00:00:00 2001 From: Mike Lin Date: Tue, 25 Dec 2012 00:21:32 -0800 Subject: [PATCH] Change over to gsl-ocaml --- lib/CamlPaml/Fit.ml | 38 ++++++++-------- lib/CamlPaml/P.ml | 6 +-- lib/CamlPaml/P.mli | 2 +- lib/CamlPaml/PhyloEM.ml | 6 +-- lib/CamlPaml/PhyloLik.ml | 60 ++++++++++++------------- lib/CamlPaml/PhyloModel.ml | 2 +- lib/CamlPaml/Q.ml | 108 ++++++++++++++++++++++----------------------- lib/CamlPaml/test.ml | 4 +- src/OmegaModel.ml | 2 +- src/PhyloCSF.ml | 2 +- src/PhyloCSFTest.ml | 2 +- 11 files changed, 116 insertions(+), 116 deletions(-) diff --git a/lib/CamlPaml/Fit.ml b/lib/CamlPaml/Fit.ml index 0a7b1f4..86a0271 100644 --- a/lib/CamlPaml/Fit.ml +++ b/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) diff --git a/lib/CamlPaml/P.ml b/lib/CamlPaml/P.ml index ebed0d1..bdc839d 100644 --- a/lib/CamlPaml/P.ml +++ b/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}) diff --git a/lib/CamlPaml/P.mli b/lib/CamlPaml/P.mli index 9bb61cc..1081d80 100644 --- a/lib/CamlPaml/P.mli +++ b/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*) diff --git a/lib/CamlPaml/PhyloEM.ml b/lib/CamlPaml/PhyloEM.ml index 8125735..081caaf 100644 --- a/lib/CamlPaml/PhyloEM.ml +++ b/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 diff --git a/lib/CamlPaml/PhyloLik.ml b/lib/CamlPaml/PhyloLik.ml index 76f6d15..6157f59 100644 --- a/lib/CamlPaml/PhyloLik.ml +++ b/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 diff --git a/lib/CamlPaml/PhyloModel.ml b/lib/CamlPaml/PhyloModel.ml index 1050b21..dc0e4bf 100644 --- a/lib/CamlPaml/PhyloModel.ml +++ b/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 diff --git a/lib/CamlPaml/Q.ml b/lib/CamlPaml/Q.ml index c0ec350..78c95bd 100644 --- a/lib/CamlPaml/Q.ml +++ b/lib/CamlPaml/Q.ml @@ -24,29 +24,29 @@ let real_of_complex ?(tol=1e-6) z = failwith (sprintf "CamlPaml.Q.real_of_complex %g+%gi" z.Complex.re z.Complex.im) z.Complex.re -let m_of_cm cm = Gsl_matrix.of_arrays (Array.map (Array.map real_of_complex) (Gsl_matrix_complex.to_arrays cm)) +let m_of_cm cm = Gsl.Matrix.of_arrays (Array.map (Array.map real_of_complex) (Gsl.Matrix_complex.to_arrays cm)) -let cm_of_m m = Gsl_matrix_complex.of_arrays (Array.map (Array.map (fun re -> { Complex.re = re; im = 0. })) (Gsl_matrix.to_arrays m)) +let cm_of_m m = Gsl.Matrix_complex.of_arrays (Array.map (Array.map (fun re -> { Complex.re = re; im = 0. })) (Gsl.Matrix.to_arrays m)) (* stupid wrappers *) let zgemm ?c a b = - let m,n = Gsl_matrix_complex.dims a - let n',p = Gsl_matrix_complex.dims b + let m,n = Gsl.Matrix_complex.dims a + let n',p = Gsl.Matrix_complex.dims b if n <> n' then invalid_arg "CamlPaml.Q.zgemm: incompatible dimensions" let c = match c with | Some c -> c - | None -> Gsl_matrix_complex.create m p - Gsl_blas.Complex.gemm ~ta:Gsl_blas.NoTrans ~tb:Gsl_blas.NoTrans ~alpha:Complex.one ~a:a ~b:b ~beta:Complex.zero ~c:c + | None -> Gsl.Matrix_complex.create m p + Gsl.Blas.Complex.gemm ~ta:Gsl.Blas.NoTrans ~tb:Gsl.Blas.NoTrans ~alpha:Complex.one ~a:a ~b:b ~beta:Complex.zero ~c:c c let zdiagm ?c d a = - let m = Gsl_vector_complex.length d - let (n,p) = Gsl_matrix_complex.dims a + let m = Gsl.Vector_complex.length d + let (n,p) = Gsl.Matrix_complex.dims a if m <> n then invalid_arg "CamlPaml.Q.diagm: incompatible dimensions" let c = match c with | Some c -> c - | None -> Gsl_matrix_complex.create m p + | None -> Gsl.Matrix_complex.create m p for i = 0 to m-1 do let d_i = d.{i} @@ -56,13 +56,13 @@ let zdiagm ?c d a = c let zinvm m = - let n,n' = Gsl_matrix_complex.dims m + let n,n' = Gsl.Matrix_complex.dims m if n <> n' then invalid_arg "CamlPaml.Q.zinvm: non-square matrix" - let p = Gsl_permut.make n - let lu = Gsl_vectmat.cmat_convert (`CM (Gsl_matrix_complex.copy m)) - ignore (Gsl_linalg.complex_LU_decomp lu p) - let m' = Gsl_vectmat.cmat_convert (`CM (Gsl_matrix_complex.create n n)) - Gsl_linalg.complex_LU_invert lu p m' + let p = Gsl.Permut.make n + let lu = Gsl.Vectmat.cmat_convert (`CM (Gsl.Matrix_complex.copy m)) + ignore (Gsl.Linalg.complex_LU_decomp lu p) + let m' = Gsl.Vectmat.cmat_convert (`CM (Gsl.Matrix_complex.create n n)) + Gsl.Linalg.complex_LU_invert lu p m' match m' with | `CM x -> x | _ -> assert false @@ -70,35 +70,35 @@ let zinvm m = module Diag = struct (* diagonalized Q = S*L*S' *) type t = { - q : Gsl_matrix.matrix; (* Q *) - s : Gsl_matrix_complex.matrix; (* S = right eigenvectors (in the columns) *) - s' : Gsl_matrix_complex.matrix; (* S' = left eigenvectors (in the rows) *) - l : Gsl_vector_complex.vector; (* diag(L) = eigenvalues *) - pi : Gsl_vector.vector; + q : Gsl.Matrix.matrix; (* Q *) + s : Gsl.Matrix_complex.matrix; (* S = right eigenvectors (in the columns) *) + s' : Gsl.Matrix_complex.matrix; (* S' = left eigenvectors (in the rows) *) + l : Gsl.Vector_complex.vector; (* diag(L) = eigenvalues *) + pi : Gsl.Vector.vector; mutable have_pi : bool; - mutable memoized_to_Pt : (float -> Gsl_matrix.matrix) option; + mutable memoized_to_Pt : (float -> Gsl.Matrix.matrix) option; mutable tol : float; } let dim q = - let (n,n') = Gsl_matrix.dims q.q + let (n,n') = Gsl.Matrix.dims q.q assert (n = n') n let of_Q ?(tol=1e-6) qm = - let qm = Gsl_matrix.of_arrays qm - let l, s = Gsl_eigen.nonsymmv ~protect:true (`M qm) + let qm = Gsl.Matrix.of_arrays qm + let l, s = Gsl.Eigen.nonsymmv ~protect:true (`M qm) let s' = zinvm s { q = qm; s = s; s' = s'; l = l; - pi = Gsl_vector.create (fst (Gsl_matrix.dims qm)); have_pi = false; memoized_to_Pt = None; tol = tol } + pi = Gsl.Vector.create (fst (Gsl.Matrix.dims qm)); have_pi = false; memoized_to_Pt = None; tol = tol } - let to_Q q = Gsl_matrix.to_arrays q.q + let to_Q q = Gsl.Matrix.to_arrays q.q let reversible q = let rev = ref true - let n = Gsl_vector_complex.length q.l + let n = Gsl.Vector_complex.length q.l for i = 0 to n-1 do if not (check_real ~tol:q.tol q.l.{i}) then rev := false !rev @@ -107,7 +107,7 @@ module Diag = struct if not q.have_pi then let min_L = ref infinity let min_Lp = ref (-1) - let n = Gsl_vector_complex.length q.l + let n = Gsl.Vector_complex.length q.l for i = 0 to n-1 do let nm_lz = Complex.norm q.l.{i} if nm_lz < !min_L then @@ -116,50 +116,50 @@ module Diag = struct assert (!min_Lp >= 0) if (abs_float !min_L) > q.tol then failwith "CamlPaml.Q.equilibrium: smallest-magnitude eigenvalue is unacceptably large; is this a valid rate matrix?" - let lev = Gsl_matrix_complex.row q.s' !min_Lp + let lev = Gsl.Matrix_complex.row q.s' !min_Lp let mass = ref Complex.zero for i = 0 to n-1 do mass := Complex.add !mass lev.{i} for i = 0 to n-1 do q.pi.{i} <- real_of_complex (Complex.div lev.{i} !mass) q.have_pi <- true - Gsl_vector.to_array q.pi + Gsl.Vector.to_array q.pi (** normalize to unity mean rate of replacement let normalize ?tol q = - let n = Gsl_vector_complex.length q.l + let n = Gsl.Vector_complex.length q.l let pi = equilibrium ?tol q let tot = ref 0. for i = 0 to n-1 do tot := !tot +. (-1.) *. pi.(i) *. q.q.{i,i} - let nq = Gsl_matrix.copy q.q - Gsl_matrix.scale q.q (1. /. !tot) - let nl = Gsl_vector_complex.copy q.l + let nq = Gsl.Matrix.copy q.q + Gsl.Matrix.scale q.q (1. /. !tot) + let nl = Gsl.Vector_complex.copy q.l for i = 0 to n-1 do nl.{i} <- Complex.div nl.{i} { Complex.re = !tot; im = 0. } { q with q = nq; l = nl } *) let scale q x = if x <= 0. then invalid_arg "CamlPaml.Q.scale: nonpositive scale factor" - let xq = Gsl_matrix.copy q.q - Gsl_matrix.scale xq x - let xl = Gsl_vector_complex.copy q.l + let xq = Gsl.Matrix.copy q.q + Gsl.Matrix.scale xq x + let xl = Gsl.Vector_complex.copy q.l let cx = { Complex.re = x; im = 0. } - for i = 0 to (Gsl_vector_complex.length xl) - 1 do + for i = 0 to (Gsl.Vector_complex.length xl) - 1 do xl.{i} <- Complex.mul xl.{i} cx { q with q = xq; l = xl; memoized_to_Pt = None } let real_to_Pt q t = if t < 0. then invalid_arg "CamlPaml.Q.to_Pt" let ct = { Complex.re = t; im = 0. } - let expLt = Gsl_vector_complex.copy q.l + let expLt = Gsl.Vector_complex.copy q.l - for i = 0 to Gsl_vector_complex.length expLt - 1 do + for i = 0 to Gsl.Vector_complex.length expLt - 1 do expLt.{i} <- Complex.exp (Complex.mul ct expLt.{i}) let sm = m_of_cm (zgemm q.s (zdiagm expLt q.s')) - let n,_ = Gsl_matrix.dims sm + let n,_ = Gsl.Matrix.dims sm for i = 0 to n-1 do let tot = ref 0. let smii = ref 1. @@ -194,7 +194,7 @@ module Diag = struct let to_Pt_gc q = q.memoized_to_Pt <- None let dPt_dt ~q ~t = - let n,_ = Gsl_matrix.dims q.q + let n,_ = Gsl.Matrix.dims q.q let (( * ),(+),(/)) = Complex.mul, Complex.add, Complex.div let exp= Complex.exp let s, l, s' = q.s, q.l, q.s' @@ -210,10 +210,10 @@ module Diag = struct (* TODO in richly parameterized models, dQ_dx is likely to be sparse. *) let dPt_dQ_dx ~q ~t ~dQ_dx = - let n,_ = Gsl_matrix.dims q.q - let dQt_dx = cm_of_m (Gsl_matrix.of_arrays dQ_dx) + let n,_ = Gsl.Matrix.dims q.q + let dQt_dx = cm_of_m (Gsl.Matrix.of_arrays dQ_dx) - if Gsl_matrix_complex.dims dQt_dx <> Gsl_matrix.dims q.q then + if Gsl.Matrix_complex.dims dQt_dx <> Gsl.Matrix.dims q.q then invalid_arg "CamlPaml.Q.dP_dx: wrong dimension of dQ_dx" let ct = { Complex.re = t; Complex.im = 0. } @@ -221,9 +221,9 @@ module Diag = struct let (( * ),(-),(/)) = Complex.mul, Complex.sub, Complex.div (* combos 1,3 or 2 (alone) -- 1,3 matches P&S? *) - Gsl_matrix_complex.scale dQt_dx ct (* 1 *) + Gsl.Matrix_complex.scale dQt_dx ct (* 1 *) - let f = Gsl_matrix_complex.create n n + let f = Gsl.Matrix_complex.create n n for i = 0 to pred n do for j = 0 to pred n do if q.l.{i} = q.l.{j} then @@ -232,13 +232,13 @@ module Diag = struct f.{i,j} <- (exp (q.l.{i} * ct) - exp (q.l.{j} * ct)) / ((q.l.{i} - q.l.{j}) * ct (* 3 *)) (* ehh not being too gentle with the allocator/GC here *) - Gsl_matrix_complex.mul_elements f (zgemm q.s' (zgemm dQt_dx q.s)) - Gsl_matrix.to_arrays (m_of_cm (zgemm q.s (zgemm f q.s'))) + Gsl.Matrix_complex.mul_elements f (zgemm q.s' (zgemm dQt_dx q.s)) + Gsl.Matrix.to_arrays (m_of_cm (zgemm q.s (zgemm f q.s'))) -let logm (m:Gsl_matrix.matrix) = - let n,n' = Gsl_matrix.dims m +let logm (m:Gsl.Matrix.matrix) = + let n,n' = Gsl.Matrix.dims m if n <> n' then invalid_arg "CamlPaml.Q.logm: non-square matrix" - let l, s = Gsl_eigen.nonsymmv ~protect:true (`M m) + let l, s = Gsl.Eigen.nonsymmv ~protect:true (`M m) let s' = zinvm s for i = 0 to n-1 do l.{i} <- Complex.log l.{i} @@ -247,7 +247,7 @@ let logm (m:Gsl_matrix.matrix) = (* Correction of a square matrix with zero row-sums but potentially negative off-diagonal entries into a valid rate matrix. As suggested by Israel, Rosenthal & Wei (2001) *) let irw2001 rawq = - let n,n' = Gsl_matrix.dims rawq + let n,n' = Gsl.Matrix.dims rawq assert (n = n') for i = 0 to n-1 do let g_i = ref (abs_float rawq.{i,i}) @@ -268,4 +268,4 @@ let irw2001 rawq = let of_P ?tol m = P.validate ?tol m - Gsl_matrix.to_arrays (irw2001 (logm m)) + Gsl.Matrix.to_arrays (irw2001 (logm m)) diff --git a/lib/CamlPaml/test.ml b/lib/CamlPaml/test.ml index d8f4e5e..2c26827 100644 --- a/lib/CamlPaml/test.ml +++ b/lib/CamlPaml/test.ml @@ -9,8 +9,8 @@ module TestInfer = struct let nt = NewickParser.parse NewickLexer.token (Lexing.from_string "(A,(B,C))") let t = T.of_newick nt - let sm0 = Gsl_matrix.of_arrays [| [| 0.8; 0.2; |]; [| 0.25; 0.75 |] |] - let sm1 = Gsl_matrix.of_arrays [| [| 0.9; 0.1; |]; [| 0.85; 0.15 |] |] + let sm0 = Gsl.Matrix.of_arrays [| [| 0.8; 0.2; |]; [| 0.25; 0.75 |] |] + let sm1 = Gsl.Matrix.of_arrays [| [| 0.9; 0.1; |]; [| 0.85; 0.15 |] |] let sms = [| sm0; sm1; sm1; sm1 |] let prior = [| 0.6; 0.4 |] diff --git a/src/OmegaModel.ml b/src/OmegaModel.ml index 217bf96..395e19c 100644 --- a/src/OmegaModel.ml +++ b/src/OmegaModel.ml @@ -154,7 +154,7 @@ let half_cauchy_lpdf ?(mode=0.0) ~scale x = log numer -. log denom let lpr_rho = half_cauchy_lpdf ~mode:1.0 ~scale:0.5 (* rho = tree_scale (as in manuscript) *) -let lpr_kappa k = log (Gsl_randist.gamma_pdf ~a:7.0 ~b:0.25 (k -. 1.0 +. epsilon_float)) +let lpr_kappa k = log (Gsl.Randist.gamma_pdf ~a:7.0 ~b:0.25 (k -. 1.0 +. epsilon_float)) (* find MAP estimates of kappa & rho *) let kr_map leaves inst = diff --git a/src/PhyloCSF.ml b/src/PhyloCSF.ml index 574a6dd..c85826c 100644 --- a/src/PhyloCSF.ml +++ b/src/PhyloCSF.ml @@ -5,7 +5,7 @@ open OptParse open Printf open CamlPaml -Gsl_error.init () +Gsl.Error.init () module SMap = Map.StringMap module SSet = Set.StringSet diff --git a/src/PhyloCSFTest.ml b/src/PhyloCSFTest.ml index 782cd15..8a916fe 100644 --- a/src/PhyloCSFTest.ml +++ b/src/PhyloCSFTest.ml @@ -7,7 +7,7 @@ open OptParse open Printf open CamlPaml -Gsl_error.init () +Gsl.Error.init () Random.self_init () module Codon = CamlPaml.Code.Codon64