From fd3a0fb9e1ea4c62072b57977eb671096563dbff Mon Sep 17 00:00:00 2001 From: Mike Lin Date: Sun, 18 Aug 2013 22:41:43 -0700 Subject: [PATCH] CamlPaml.PhyloModel.P14n.update: simplify 'inheritance' of root prior distribution for reversible models --- lib/CamlPaml/PhyloModel.ml | 34 ++++++++++++++++++---------------- lib/CamlPaml/PhyloModel.mli | 15 +++++++++++++-- 2 files changed, 31 insertions(+), 18 deletions(-) diff --git a/lib/CamlPaml/PhyloModel.ml b/lib/CamlPaml/PhyloModel.ml index dc0e4bf..30f3e55 100644 --- a/lib/CamlPaml/PhyloModel.ml +++ b/lib/CamlPaml/PhyloModel.ml @@ -6,7 +6,7 @@ type t = { tree : T.t; qms : Q.Diag.t array; pms : P.matrix array; - prior : float array + prior : float array option } let make ?prior t qms = @@ -16,25 +16,28 @@ let make ?prior t qms = if Q.Diag.dim qms.(i) <> Q.Diag.dim qms.(0) || T.branch t i < 0. then invalid_arg "CamlPaml.PhyloModel.make" let pms = Array.init (T.size t - 1) (fun br -> Q.Diag.to_Pt qms.(br) (T.branch t br)) let prior = match prior with - | Some pr -> Array.copy pr + | Some pr -> + if Array.length pr <> Q.Diag.dim qms.(0) then invalid_arg "CamlPaml.PhyloModel.make" + Some (Array.copy pr) | None -> - let q = qms.(snd (T.children t (T.root t))) - if not (Q.Diag.reversible q) then invalid_arg "CamlPaml.PhyloModel.make" - Q.Diag.equilibrium q - { tree = T.copy t; qms = qms; pms = pms; prior = prior } + if not (Q.Diag.reversible qms.(snd (T.children t (T.root t)))) then invalid_arg "CamlPaml.PhyloModel.make" + None + { tree = T.copy t; qms; pms; prior } -let tree { tree = t } = T.copy t -let q { qms = qms } br = qms.(br) -let p { pms = pms } br = pms.(br) (* Array.map Array.copy pms.(br) *) -let prior { prior = rp } = Array.copy rp +let tree { tree } = T.copy tree +let q { qms } br = qms.(br) +let p { pms } br = pms.(br) (* Array.map Array.copy pms.(br) *) +let prior { tree; qms; prior } = match prior with + | Some pr -> Array.copy pr + | None -> Q.Diag.equilibrium qms.(snd (T.children tree (T.root tree))) (* q was verified reversible in [make], above*) -let prepare_lik ?workspace m leaves = PhyloLik.prepare ?workspace m.tree m.pms m.prior leaves +let prepare_lik ?workspace m leaves = PhyloLik.prepare ?workspace m.tree m.pms (prior m) leaves 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 root_chooser = Tools.random_chooser ~checksum:checksum m.prior + let root_chooser = Tools.random_chooser ~checksum:checksum (prior m) fun ?root ?a () -> let t = m.tree let a = match a with @@ -138,10 +141,9 @@ module P14n = struct pms_changed := true 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 - | None -> inst.model.prior - let model = if !pms_changed then make ~prior newtree newq else { inst.model with prior = prior } + let model = + if !pms_changed then make ?prior newtree newq + else { inst.model with prior = prior } { inst with model = model; q_settings = (match q_settings with Some set -> Array.copy set | None -> inst.q_settings); tree_settings = (match tree_settings with Some set -> Array.copy set | None -> inst.tree_settings) } diff --git a/lib/CamlPaml/PhyloModel.mli b/lib/CamlPaml/PhyloModel.mli index 0bed36f..7c184a3 100644 --- a/lib/CamlPaml/PhyloModel.mli +++ b/lib/CamlPaml/PhyloModel.mli @@ -61,7 +61,7 @@ module P14n : sig (** in an instance of a p14n we have settings for the variables, thus determining a fully parameterized model *) type instance - (** instantiate the model by giving settings for all the variables in its p14n. *) + (** instantiate the model by giving settings for all the variables in its p14n. If prior is not provided, it is determined as in [PhyloModel.make], above. *) val instantiate : ?prior:(float array) -> model_p14n -> q_settings:(float array) -> tree_settings:(float array) -> instance val model : instance -> t @@ -69,7 +69,18 @@ module P14n : sig val q_settings : instance -> float array val tree_settings : instance -> float array - (** return a copy of the instance with a subset of the settings replaced (possibly reusing computed information in the original instance that is not affected by the changed parameters) *) + (** return a copy of the instance with a subset of the settings replaced (possibly reusing computed information in the original instance that is not affected by the changed parameters) + + When [prior] is not provided, the prior distribution for the updated model + is determined as follows. If the prior of the provided instance (or any of + its 'ancestral' instances) had been explicitly set in [instantiate], the + updated model will have the same previously set prior. If the provided + instance's prior was instead determined based on rate matrix equilibrium + frequencies (i.e. no prior was provided to [instantiate]), the new + instance's prior will be determined based on the rate matrix equilibrium + frequencies in the new model; thus, if [q_settings] changes as a result of + the update, the prior may also change. + *) val update : ?prior:(float array) -> ?q_settings:(float array) -> ?tree_settings:(float array) -> instance -> instance