From bafc397a61b83f768a009af057a58ab9e5e937d6 Mon Sep 17 00:00:00 2001 From: Michael Lin Date: Mon, 7 Feb 2011 17:36:03 -0500 Subject: [PATCH] CamlPaml: PhyloModel.P14n checks variable domains; correction to Fit.maximizer documentation --- lib/CamlPaml/Fit.ml | 22 +++++++++++++--------- lib/CamlPaml/Fit.mli | 5 ++++- lib/CamlPaml/PhyloModel.ml | 20 ++++++++++++-------- ocaml+twt | 1 - 4 files changed, 29 insertions(+), 19 deletions(-) delete mode 160000 ocaml+twt diff --git a/lib/CamlPaml/Fit.ml b/lib/CamlPaml/Fit.ml index e8042c1..8b803da 100644 --- a/lib/CamlPaml/Fit.ml +++ b/lib/CamlPaml/Fit.ml @@ -96,21 +96,25 @@ class multi_maximizer f df init = let make_multi_maximizer ~f ~df ~init = (new multi_maximizer f df init) type domain = Real | Pos | Neg | NonPos | NonNeg | Probability - + +let check_domain domain x = + match domain with + | Real when x=x -> true + | Pos when x > 0. -> true + | Neg when x < 0. -> true + | NonPos when x <= 0. -> true + | NonNeg when x >= 0. -> true + | Probability when x >= 0. && x <= 1. -> true + | _ -> false + + exception False let enforce_domain ?(rejection_value=neg_infinity) f domains = fun x -> if (Array.length x) <> Array.length domains then invalid_arg "CamlPaml.Fit.enforce_domains: length mismatch" try Array.iteri - fun i domain -> - match domain with - | Pos when x.(i) <= 0. -> raise False - | Neg when x.(i) >= 0. -> raise False - | NonPos when x.(i) > 0. -> raise False - | NonNeg when x.(i) < 0. -> raise False - | Probability when x.(i) < 0. || x.(i) > 1. -> raise False - | _ -> () + fun i domain -> if not (check_domain domain x.(i)) then raise False domains f x with diff --git a/lib/CamlPaml/Fit.mli b/lib/CamlPaml/Fit.mli index e4e5bc7..ad7b2a0 100644 --- a/lib/CamlPaml/Fit.mli +++ b/lib/CamlPaml/Fit.mli @@ -4,7 +4,7 @@ (** a [maximizer] iteratively narrows down a local maximum of a one-dimensional function [f(x)] *) class type maximizer = object - (** current estimate of the maximum [max_x f(x)] *) + (** current estimate of the locally maximizing setting, [argmax_x f(x)] *) method maximum : unit -> float (** current lower and upper bounds on [x] around the local maximum *) @@ -47,5 +47,8 @@ val make_multi_maximizer : f:(float array -> float) -> df:(float array -> float type domain = Real | Pos | Neg | NonPos | NonNeg | Probability +(** [check_domain domain x] returns [true] if [x] is in [domain], false otherwise *) +val check_domain : domain -> float -> bool + (** wrap the function [f(x)] to prevent the [multi_maximizer] from searching outside of the specified domain for each element of the vector input. *) val enforce_domain : ?rejection_value:float -> (float array -> float) -> (domain array) -> (float array -> float) diff --git a/lib/CamlPaml/PhyloModel.ml b/lib/CamlPaml/PhyloModel.ml index 481b63d..25b47b3 100644 --- a/lib/CamlPaml/PhyloModel.ml +++ b/lib/CamlPaml/PhyloModel.ml @@ -77,13 +77,17 @@ module P14n = struct tot := Expr.Add (q.(i).(j), !tot) q.(i).(i) <- Expr.simplify (Expr.Sub (Expr.Val 0., !tot)) - let instantiate_tree shape exprs settings = + let instantiate_tree shape exprs domains settings = + Array.iteri (fun i domain -> if not (Fit.check_domain domain settings.(i)) then invalid_arg ("CamlPaml.P14n.instantiate_tree: domain violation on variable " ^ (string_of_int i))) domains + let tree = T.copy shape for br = 0 to T.root tree - 1 do T.put_branch tree br (Expr.eval exprs.(br) settings) tree - let instantiate_q exprs scale_expr settings = + let instantiate_q exprs scale_expr domains settings = + Array.iteri (fun i domain -> if not (Fit.check_domain domain settings.(i)) then invalid_arg ("CamlPaml.P14n.instantiate_q: domain violation on variable " ^ (string_of_int i))) domains + let scale = Expr.eval scale_expr settings if scale <= 0. then @@ -95,7 +99,7 @@ module P14n = struct q - let instantiate_qs p14ns scale_p14ns settings = + let instantiate_qs p14ns scale_p14ns domains settings = let previous = ref [] (* memoized results from previous branches...I'm assuming the number of independent rate matrix parameterizations to be sublinear in the size of the tree, otherwise this memoization is quadratic... *) Array.init (Array.length p14ns) fun br -> @@ -109,13 +113,13 @@ module P14n = struct q with | Not_found -> - let q = instantiate_q p14ns.(br) scale_p14ns.(br) settings + let q = instantiate_q p14ns.(br) scale_p14ns.(br) domains settings previous := (p14ns.(br),scale_p14ns.(br),q) :: !previous q let instantiate ?prior p14n ~q_settings ~tree_settings = - let qms = instantiate_qs p14n.q_p14ns p14n.q_scale_p14ns q_settings - let tree = instantiate_tree p14n.tree_shape p14n.tree_p14n tree_settings + let qms = instantiate_qs p14n.q_p14ns p14n.q_scale_p14ns p14n.q_domains q_settings + let tree = instantiate_tree p14n.tree_shape p14n.tree_p14n p14n.tree_domains tree_settings { model = make ?prior tree qms; p14n = p14n; q_settings = Array.copy q_settings; tree_settings = Array.copy tree_settings } let update ?prior ?q_settings ?tree_settings inst = @@ -123,12 +127,12 @@ module P14n = struct let newq = match q_settings with | Some q_settings -> pms_changed := true - instantiate_qs inst.p14n.q_p14ns inst.p14n.q_scale_p14ns q_settings + instantiate_qs inst.p14n.q_p14ns inst.p14n.q_scale_p14ns inst.p14n.q_domains q_settings | None -> inst.model.qms let newtree = match tree_settings with | Some tree_settings -> pms_changed := true - instantiate_tree inst.p14n.tree_shape inst.p14n.tree_p14n tree_settings + instantiate_tree inst.p14n.tree_shape inst.p14n.tree_p14n inst.p14n.tree_domains tree_settings | None -> inst.model.tree let prior = match prior with | Some ar -> Array.copy ar diff --git a/ocaml+twt b/ocaml+twt deleted file mode 160000 index 052ee1f..0000000 --- a/ocaml+twt +++ /dev/null @@ -1 +0,0 @@ -Subproject commit 052ee1fc64866e9ae380cbd28df16b07fcf83e82