Commit 5ad509bf authored by Louis Duchemin's avatar Louis Duchemin
Browse files

Mutsel simulator with CpG entrypoint in Codepi

Also matching Phylogenetics RNG implementation using Gsl
parent 9fee54e9
open Core_kernel
open Bistro
open Codepitk
let simulation ?branch_scale ?seed ~n_h0 ~n_ha ~ne_s:(ne_s0, ne_s1) ~gBGC:(gBGC0, gBGC1) ~tree ~fitness_profiles () =
let simulation ?branch_scale ~rng ~n_h0 ~n_ha
~ne_s:(ne_s0, ne_s1) ~gBGC:(gBGC0, gBGC1)
?rate_CpG ~tree ~fitness_profiles () =
let f = fun%workflow () ->
let open Phylogenetics in
let () = Option.iter ~f:Random.init [%param seed] in
let n_h0 = [%param n_h0] in
let n_ha = [%param n_ha] in
let ne_s0 = [%param ne_s0] in
......@@ -12,33 +15,42 @@ let simulation ?branch_scale ?seed ~n_h0 ~n_ha ~ne_s:(ne_s0, ne_s1) ~gBGC:(gBGC0
let gBGC0 = [%param gBGC0] in
let gBGC1 = [%param gBGC1] in
let branch_scale = [%param branch_scale] in
let rate_CpG = [%param rate_CpG] in
let tree =
Codepitk.Convergence_tree.from_file [%path tree]
Convergence_tree.from_file [%path tree]
|> Rresult.R.get_ok
in
let fitness_profiles =
Codepitk.Profile_tsv.(
Profile_tsv.(
read [%path fitness_profiles]
|> to_fitness
|> Array.map ~f: Amino_acid.Vector.of_array_exn
) in
let tree = Option.value_map branch_scale ~default:tree
~f:(fun scale -> Codepitk.Mutsel_simulator.tree_rescale_branch_length tree ~scale)
~f:(fun scale -> Mutsel_simulation.tree_rescale_branch_length tree ~scale)
in
Codepitk.Mutsel_simulator.Site_independent.make
?seed ~n_h0 ~n_ha ~gBGC:(gBGC0, gBGC1)
~ne_s:(ne_s0, ne_s1) ~tree ~fitness_profiles ()
match rate_CpG with
|Some rate_CpG ->
Mutsel_simulation.make_with_CpG_hypermutability
~rng ~n_h0 ~n_ha ~ne_s:(ne_s0, ne_s1) ~gBGC:(gBGC0, gBGC1)
~rate_CpG ~tree ~fitness_profiles ()
|None ->
Mutsel_simulation.make_site_independent
~rng ~n_h0 ~n_ha ~gBGC:(gBGC0, gBGC1)
~ne_s:(ne_s0, ne_s1) ~tree ~fitness_profiles ()
in
Workflow.plugin ~descr:"simulator.simulation" f
let alignment_of_simulation sim =
let f = fun%workflow dest ->
let sim : Codepitk.Mutsel_simulator.Site_independent.t =
[%eval sim]
let sim : Mutsel_simulation.t = [%eval sim] in
let tree, sequences = match sim with
| Site_independent sim -> sim.tree, sim.sequences
| With_CpG_hypermutability sim -> sim.tree, sim.sequences
in
let species_name = Phylogenetics.Tree.leaves sim.tree in
let species_name = Phylogenetics.Tree.leaves tree in
Out_channel.with_file dest ~f:(fun oc ->
List.iter2_exn species_name sim.sequences ~f:(fun description sequence ->
List.iter2_exn species_name sequences ~f:(fun description sequence ->
fprintf oc ">%s\n%s\n" description (sequence :> string)
)
)
......@@ -48,7 +60,7 @@ let alignment_of_simulation sim =
let fitness_histogram sim =
let f = fun%workflow dest ->
let open Phylogenetics in
let sim : Codepitk.Mutsel_simulator.Site_independent.t =
let sim : Mutsel_simulation.Site_independent.t =
[%eval sim]
in
let params = Array.append sim.h0_params sim.ha_params in
......
......@@ -4,20 +4,21 @@ open Codepitk
val simulation :
?branch_scale:float ->
?seed:int ->
rng:Gsl.Rng.t ->
n_h0:int ->
n_ha:int ->
ne_s:float * float ->
gBGC:float * float ->
?rate_CpG: float ->
tree:nhx file ->
fitness_profiles:#text file ->
unit ->
Mutsel_simulator.Site_independent.t workflow
Mutsel_simulation.t workflow
val alignment_of_simulation :
Mutsel_simulator.Site_independent.t workflow ->
Mutsel_simulation.t workflow ->
nucleotide_fasta file
val fitness_histogram :
Mutsel_simulator.Site_independent.t workflow ->
Mutsel_simulation.Site_independent.t workflow ->
pdf file
......@@ -74,9 +74,10 @@ module Mutsel_query = struct
ne_s : float * float ;
gBGC : float * float ;
seed : int ;
rate_CpG : float option;
}
let make ?(branch_scale = 1.) ?(ne_s = 1., 1.) ?(gBGC = 0., 0.) ?(seed = 0) ~tree ~profiles ~n_h0 ~n_ha () = {
let make ?(branch_scale = 1.) ?(ne_s = 1., 1.) ?(gBGC = 0., 0.) ?rate_CpG ?(seed = 0) ~tree ~profiles ~n_h0 ~n_ha () = {
tree ;
profiles ;
n_h0 ;
......@@ -85,12 +86,15 @@ module Mutsel_query = struct
gBGC ;
branch_scale ;
seed : int ;
rate_CpG;
}
let simulation { n_h0 ; n_ha ; profiles ; ne_s ; gBGC ; branch_scale ; seed ; tree ; _ } =
let simulation { n_h0 ; n_ha ; profiles ; ne_s ; gBGC ; branch_scale ; seed ; tree ; rate_CpG ; _ } =
let rng = Gsl.Rng.(make (default ())) in
Gsl.Rng.set rng (Nativeint.of_int seed);
let tree = tree_workflow tree in
let fitness_profiles = Workflow.input profiles in
Mutsel_simulator.simulation ~branch_scale ~n_ha ~n_h0 ~ne_s ~gBGC ~tree ~seed ~fitness_profiles ()
Mutsel_simulator.simulation ~branch_scale ~n_ha ~n_h0 ~ne_s ~gBGC ?rate_CpG ~tree ~rng ~fitness_profiles ()
let nucleotide_alignment p =
simulation p
......@@ -116,18 +120,31 @@ module Mutsel = struct
method_labels : string list ;
method_outputs : float option array list ;
average_precision : (float * (float * float)) list ;
profiles : (float array * float array) array ;
profiles : (float array * float array) array option;
ancestral_counts : int Amino_acid.table array ;
convergent_counts : int Amino_acid.table array ;
}
let counts seqs i =
let module Codon = Codon.Universal_genetic_code.NS in
Amino_acid.Table.init (fun aa ->
let aa = Amino_acid.to_char aa in
List.count seqs ~f:(fun s ->
let codon_str = String.sub (s : Dna.t :> string) ~pos:(i * 3) ~len:3 in
let codon = match Codon.of_string codon_str with
| Some c -> c
| None -> assert false
in
Char.equal (Amino_acid.to_char (Codon.aa_of_codon codon)) aa)
)
let benchmark q methods =
let f = fun%workflow () ->
let open Phylogenetics in
let open Codepitk in
let open Codepitk.Mutsel_simulator.Site_independent in
(* let open Codepitk.Mutsel_simulator.Site_independent in *)
let module Codon = Codon.Universal_genetic_code.NS in
let sim : t = [%eval simulation q] in
let sim = [%eval simulation q] in
let result_paths = [%eval Bistro.Workflow.path_list (List.map methods ~f:(fun f -> f q))] in
let results =
List.map result_paths ~f:Cpt.of_file
......@@ -136,8 +153,7 @@ module Mutsel = struct
|> List.concat_map ~f:Cpt.columns
in
let method_labels, method_outputs = List.unzip results in
let n_h0 = Array.length sim.h0_params in
let n_ha = Array.length sim.ha_params in
let n_h0, n_ha = q.n_h0, q.n_ha in
let nsites = n_h0 + n_ha in
let amino_acid_vector_of_codon_vector xs =
Amino_acid.Vector.init (fun aa ->
......@@ -153,28 +169,20 @@ module Mutsel = struct
|> amino_acid_vector_of_codon_vector
|> Amino_acid.Vector.to_array
in
let profiles =
Array.append sim.h0_params sim.ha_params
|> Array.map ~f:(fun (p1, p2) ->
compute_profile p1, compute_profile p2
)
in
let counts seqs i =
Amino_acid.Table.init (fun aa ->
let aa = Amino_acid.to_char aa in
List.count seqs ~f:(fun s ->
let codon_str = String.sub (s : Dna.t :> string) ~pos:(i * 3) ~len:3 in
let codon = match Codon.of_string codon_str with
| Some c -> c
| None -> assert false
in
Char.equal (Amino_acid.to_char (Codon.aa_of_codon codon)) aa)
let profiles = match (sim : Codepitk.Mutsel_simulation.t) with
| Site_independent sim ->
Some (
Array.append sim.h0_params sim.ha_params
|> Array.map ~f:(fun (p1, p2) ->
compute_profile p1, compute_profile p2
)
)
| With_CpG_hypermutability _ -> None
in
let collect_counts cond =
let species = Convergence_tree.leaves sim.tree in
let species = Convergence_tree.leaves (Mutsel_simulation.tree sim) in
let seqs =
List.map2_exn sim.sequences species ~f:(fun s (_, cond_s) ->
List.map2_exn (Mutsel_simulation.sequences sim) species ~f:(fun s (_, cond_s) ->
if Poly.equal cond cond_s then Some s else None
)
|> List.filter_opt
......@@ -208,12 +216,16 @@ module Mutsel = struct
n_h0 ; n_ha ; method_outputs ; profiles ;
ancestral_counts ; convergent_counts } = [%eval b] in
let open OCamlR_base in
let collect_profiles sel =
let collect_profiles profiles sel =
Array.map profiles ~f:sel
|> Numeric.Matrix.of_arrays
in
let ancestral_profiles = collect_profiles fst in
let convergent_profiles = collect_profiles snd in
let maybe_profiles = Option.value_map profiles ~default:[] ~f:(fun profiles ->
[
Some "ancestral_profiles", Numeric.Matrix.to_sexp (collect_profiles profiles fst) ;
Some "convergent_profiles", Numeric.Matrix.to_sexp (collect_profiles profiles snd) ;
])
in
let collect_counts c =
Array.map c ~f:(fun c -> (c : int Amino_acid.table :> int array))
|> Integer.Matrix.of_arrays
......@@ -238,15 +250,14 @@ module Mutsel = struct
l, `Numeric (Numeric.of_array_opt r)
) in
let results = Dataframe.create columns in
let data = List_.create [
Some "results", Dataframe.to_sexp results ;
Some "oracle", Logical.to_sexp oracle ;
Some "ancestral_profiles", Numeric.Matrix.to_sexp ancestral_profiles ;
Some "convergent_profiles", Numeric.Matrix.to_sexp convergent_profiles ;
Some "ancestral_counts", Integer.Matrix.to_sexp ancestral_counts ;
Some "convergent_counts", Integer.Matrix.to_sexp convergent_counts ;
Some "auc_estimates", Dataframe.to_sexp auc_estimates ;
]
let data = List_.create (
[
Some "results", Dataframe.to_sexp results ;
Some "oracle", Logical.to_sexp oracle ;
Some "ancestral_counts", Integer.Matrix.to_sexp ancestral_counts ;
Some "convergent_counts", Integer.Matrix.to_sexp convergent_counts ;
Some "auc_estimates", Dataframe.to_sexp auc_estimates ;
] @ maybe_profiles)
in
saveRDS ~file:dest (List_.to_sexp data)
in
......
......@@ -29,6 +29,7 @@ module Mutsel : sig
?branch_scale:float ->
?ne_s:float * float ->
?gBGC:float * float ->
?rate_CpG: float ->
?seed:int ->
tree:tree ->
profiles:string ->
......@@ -39,7 +40,7 @@ module Mutsel : sig
val simulation :
query ->
Codepitk.Mutsel_simulator.Site_independent.t workflow
Codepitk.Mutsel_simulation.t workflow
type benchmark = {
......@@ -48,7 +49,7 @@ module Mutsel : sig
method_labels : string list ;
method_outputs : float option array list ;
average_precision : (float * (float * float)) list ;
profiles : (float array * float array) array ;
profiles : (float array * float array) array option ;
ancestral_counts : int Phylogenetics.Amino_acid.table array ;
convergent_counts : int Phylogenetics.Amino_acid.table array ;
}
......
......@@ -81,9 +81,9 @@ module Site_independent = struct
ne_s : float * float ;
}
let alignment tree ~root param =
let alignment rng tree ~root param =
List.init (Array.length root) ~f:(fun i ->
Sim.site_gillespie_direct tree ~root:root.(i) ~param:(param i)
Sim.site_gillespie_direct rng tree ~root:root.(i) ~param:(param i)
|> Tree.leaves
|> List.map ~f:Codon.Universal_genetic_code.NS.to_string
)
......@@ -91,8 +91,7 @@ module Site_independent = struct
|> List.map ~f:String.concat
|> List.map ~f:Dna.of_string_unsafe
let make ?seed ~n_h0 ~n_ha ~ne_s ~gBGC ~tree ~fitness_profiles () =
let () = Option.iter ~f:Random.init seed in
let make ~rng ~n_h0 ~n_ha ~ne_s ~gBGC ~tree ~fitness_profiles () =
let h0_profiles, ha_profiles = make_profiles n_h0 n_ha fitness_profiles in
let base_param, h0_params, ha_params =
make_params ~gBGC ~ne_s ~fitness_profiles ~h0_profiles ~ha_profiles in
......@@ -105,8 +104,8 @@ module Site_independent = struct
|> Mutsel.NSCodon.Table.of_array_exn
)
in
let root = Sim.hmm0 ~len:(n_h0 + n_ha) ~dist:(Array.get root_dists) in
let sequences = alignment tree ~root params in
let root = Sim.hmm0 rng ~len:(n_h0 + n_ha) ~dist:(Array.get root_dists) in
let sequences = alignment rng tree ~root params in
{ sequences ; fitness_profiles ; h0_profiles ; ha_profiles ;
gBGC ; ne_s ; h0_params ; ha_params ; base_param ; tree}
end
......@@ -159,4 +158,26 @@ module With_CpG_hypermutability = struct
{ sequences ; fitness_profiles ; h0_profiles ; ha_profiles ;
gBGC ; ne_s ; tree ; rate_CpG}
end
\ No newline at end of file
end
type t =
Site_independent of Site_independent.t
| With_CpG_hypermutability of With_CpG_hypermutability.t
let make_with_CpG_hypermutability ~rng ~n_h0 ~n_ha ~ne_s ~gBGC ~rate_CpG ~tree ~fitness_profiles () =
With_CpG_hypermutability (
With_CpG_hypermutability.make ~rng ~n_h0 ~n_ha ~ne_s ~gBGC ~rate_CpG ~tree ~fitness_profiles ()
)
let make_site_independent ~rng ~n_h0 ~n_ha ~ne_s ~gBGC ~tree ~fitness_profiles () =
Site_independent (
Site_independent.make ~rng ~n_h0 ~n_ha ~ne_s ~gBGC ~tree ~fitness_profiles ()
)
let tree = function
Site_independent sim -> sim.tree
| With_CpG_hypermutability sim -> sim.tree
let sequences = function
Site_independent sim -> sim.sequences
| With_CpG_hypermutability sim -> sim.sequences
\ No newline at end of file
......@@ -18,7 +18,7 @@ module Site_independent : sig
}
val make :
?seed:int ->
rng:Gsl.Rng.t ->
n_h0:int ->
n_ha:int ->
ne_s:float * float ->
......@@ -52,5 +52,36 @@ module With_CpG_hypermutability : sig
fitness_profiles:Amino_acid.vector array ->
unit ->
t
end
type t =
Site_independent of Site_independent.t
| With_CpG_hypermutability of With_CpG_hypermutability.t
val make_with_CpG_hypermutability :
rng:Gsl.Rng.t ->
n_h0:int ->
n_ha:int ->
ne_s:float * float ->
gBGC:float * float ->
rate_CpG:float ->
tree:Convergence_tree.t ->
fitness_profiles:Amino_acid.vector array ->
unit ->
t
val make_site_independent :
rng:Gsl.Rng.t ->
n_h0:int ->
n_ha:int ->
ne_s:float * float ->
gBGC:float * float ->
tree:Convergence_tree.t ->
fitness_profiles:Amino_acid.vector array ->
unit ->
t
val tree : t -> Convergence_tree.t
val sequences : t -> Dna.t list
\ No newline at end of file
......@@ -75,8 +75,8 @@ module Evolution_model = struct
let exp_matrix = Amino_acid.Matrix.(expm (scal_mul t m)) in
Matrix.robust_equal ~tol:1e-10 diag_exp_matrix (exp_matrix :> mat)
let%test "Matrix exponential through diagonalisation matches naive implementation" =
test_diagonal_matrix_exponential ()
let%test "Matrix exponential through diagonalisation matches naive implementation" =
test_diagonal_matrix_exponential ()
end
let choose_aa p =
......@@ -119,6 +119,7 @@ module type S = sig
float * param
val simulate_site :
rng : Gsl.Rng.t ->
exchangeability_matrix:Amino_acid.matrix ->
stationary_distribution:Amino_acid.vector ->
(_, _, branch_info) Tree.t ->
......@@ -171,6 +172,7 @@ module type S = sig
Model2.param * param * likelihood_ratio_test
val simulate_site :
rng : Gsl.Rng.t ->
exchangeability_matrix:Rate_matrix.Amino_acid.t ->
scale:float ->
stationary_distribution0:Amino_acid.vector ->
......@@ -241,14 +243,14 @@ module Make(Branch_info : Branch_info)(Leaf_info : Leaf_info)(Site : Site with t
let ll, vec = inner_maximum_log_likelihood ?debug ~exchangeability_matrix ~stationary_distribution tree site in
ll, decode_vec vec
let simulate_site ~exchangeability_matrix ~stationary_distribution tree ~param:scale =
let simulate_site ~rng ~exchangeability_matrix ~stationary_distribution tree ~param:scale =
let root = choose_aa stationary_distribution in
let p = {
Evolution_model.stationary_distribution ; scale ;
exchangeability_matrix ;
}
in
Simulator.site_gillespie_first_reaction tree ~root ~param:(Fn.const p)
Simulator.site_gillespie_first_reaction rng tree ~root ~param:(Fn.const p)
end
......@@ -473,7 +475,7 @@ module Make(Branch_info : Branch_info)(Leaf_info : Leaf_info)(Site : Site with t
decode_vec schema3 p3,
lrt
let simulate_site ~exchangeability_matrix ~scale ~stationary_distribution0 ~stationary_distribution1 tree =
let simulate_site ~rng ~exchangeability_matrix ~scale ~stationary_distribution0 ~stationary_distribution1 tree =
let param b =
evolution_model_param
exchangeability_matrix
......@@ -481,7 +483,7 @@ module Make(Branch_info : Branch_info)(Leaf_info : Leaf_info)(Site : Site with t
(Branch_info.condition b)
in
let root = choose_aa stationary_distribution0 in
Simulator.site_gillespie_first_reaction tree ~root ~param
Simulator.site_gillespie_first_reaction rng tree ~root ~param
end
end
......@@ -530,12 +532,13 @@ module Implementation_check = struct
open Pack
let likelihood_plot_demo (wag : Wag.t) =
let rng = Gsl.Rng.(make (default ())) in
let tree = pair_tree ~branch_length1:1. ~branch_length2:1. ~npairs:100 in
let root = choose_aa wag.freqs in
let true_scale = 1. in
let p = Evolution_model.param_of_wag wag true_scale in
let site =
Simulator.site_gillespie_first_reaction tree ~root ~param:(Fn.const p)
Simulator.site_gillespie_first_reaction rng tree ~root ~param:(Fn.const p)
|> Site.of_simulation
in
let ll, scale_hat =
......@@ -562,12 +565,14 @@ module Implementation_check = struct
OCamlR_graphics.plot ~x ~y ()
let lrt_1_vs_2_null_simulation ?(seed = 31415926535897931) ?mode ?(nb_simulations = 1_000) (wag : Wag.t) =
Owl_stats_prng.init seed ;
let rng = Gsl.Rng.(make (default ())) in
Gsl.Rng.set rng (Nativeint.of_int seed);
let tree = pair_tree ~branch_length1:1. ~branch_length2:1. ~npairs:30 in
let true_scale = 1. in
let f _ =
let simulation =
Model1.simulate_site
~rng
~exchangeability_matrix:wag.rate_matrix
~stationary_distribution:wag.freqs
~param:true_scale tree
......@@ -579,13 +584,15 @@ module Implementation_check = struct
Array.init nb_simulations ~f
let lrt_2_vs_3_null_simulation ?(seed = 31415926535897931) ?mode ?(alpha = 0.1) ?(nb_simulations = 1_000) (wag : Wag.t) =
Owl_stats_prng.init seed ;
let rng = Gsl.Rng.(make (default ())) in
Gsl.Rng.set rng (Nativeint.of_int seed);
let tree = pair_tree ~branch_length1:1. ~branch_length2:1. ~npairs:30 in
let true_scale = 1. in
let f _ =
let stationary_distribution = simulate_profile alpha in
let simulation =
Model1.simulate_site
~rng
~exchangeability_matrix:wag.rate_matrix
~stationary_distribution
~param:true_scale tree
......
......@@ -39,6 +39,7 @@ module type S = sig
float * param
val simulate_site :
rng : Gsl.Rng.t ->
exchangeability_matrix:Amino_acid.matrix ->
stationary_distribution:Amino_acid.vector ->
(_, _, branch_info) Tree.t ->
......@@ -91,6 +92,7 @@ module type S = sig
Model2.param * param * likelihood_ratio_test
val simulate_site :
rng : Gsl.Rng.t ->
exchangeability_matrix:Rate_matrix.Amino_acid.t ->
scale:float ->
stationary_distribution0:Amino_acid.vector ->
......
......@@ -28,7 +28,7 @@ let _oneline_rodent = {
ne_s = 4., 4. ;
}
let _rubisco = {
let rubisco = {
label = "rubisco" ;
tree = Rubisco_dataset.(Path "data/rubisco" |> Query.tree ~branch_length_unit:`Amino_acid) ;
rooted = false ;
......@@ -36,7 +36,7 @@ let _rubisco = {
ne_s = 4., 4. ;
}
let _orthomam_echolocation () = {
let _orthomam_echolocation = {
label = "orthomam_echolocation" ;
tree = (
Orthomam.tree_of_db
......@@ -49,20 +49,33 @@ let _orthomam_echolocation () = {
ne_s = 4., 4. ;
}
let pvals_of_cpt convergence_table pvalue_column =
let cpt = Cpt.of_file_exn convergence_table in
Cpt.get_col_exn cpt pvalue_column
|> Array.filter_opt
|> Array.map ~f:(fun p -> 1. -. p)
let hist_pval (convergence_table : cpt file) (pvalue_column:string): pdf file =
let f = fun%workflow dest ->
let cpt = Cpt.of_file_exn [%path convergence_table] in
let pvals =
Cpt.get_col_exn cpt pvalue_column
|> Array.filter_opt
|> Array.map ~f:(fun p -> 1. -. p)
in
let pvals = pvals_of_cpt [%path convergence_table] pvalue_column in
OCamlR_grDevices.pdf dest;
OCamlR_graphics.hist ~breaks:(`n 60) pvals
|> fun _ -> OCamlR_grDevices.dev_off ()
in
Bistro.Workflow.path_plugin ~descr:"Calibration.hist_pval" f
(* let quantiles_pval convergence_table pvalue_column =
let module R = Bistro_utils.R_script in
let f = fun%workflow dest ->
let pvals = pvals_of_cpt [%path convergence_table] pvalue_column
|> Array.to_list in
[%script {|
library(tidyverse)
pvals = <<< R.(make [floats pvals]) >>>
|}] in
let script = Workflow.path_plugin ~descr:"calibration.quantiles_pval_script" f in
R.workflow ~descr:"calibration.quantiles_pval" f *)
module Under_mutsel = struct
module Pipeline = Codepi.Simulation_pipeline.Mutsel
......@@ -136,8 +149,9 @@ module Under_mutsel = struct
rds$auc_estimates %>%
ggplot(aes(x=method, y=estimate)) +
geom_pointrange(aes(ymin=lower_bound, ymax=upper_bound)) +
ylab("AUC") +
theme_bw() +
ggtitle("<<< dep rds >>>")
theme(axis.text.x = element_text(angle=45, hjust=1))
dev.off()
|}] in
Bistro_utils.R_script.workflow ~descr:"calibration.under_mutsel.plot_benchmark" script
......@@ -146,16 +160,16 @@ end
module Under_tdg09 = struct
let site_simulation ?(alpha = 1.) ?(scale = 1.) (wag : Phylogenetics.Wag.t) tree =
let site_simulation ?(alpha = 1.) ?(scale = 1.) ~rng (wag : Phylogenetics.Wag.t) tree =
let open Codepitk.Tdg09.Pack in
let exchangeability_matrix = wag.rate_matrix in
let stationary_distribution = simulate_profile alpha in
let param = scale in
Model1.simulate_site ~exchangeability_matrix ~stationary_distribution ~param tree
Model1.simulate_site ~rng ~exchangeability_matrix ~stationary_distribution ~param tree
let simulation ?alpha wag tree nsites =
let simulation ~rng ?alpha wag tree nsites =
Array.init nsites ~f:(fun _ ->
site_simulation wag ?alpha tree
site_simulation ~rng wag ?alpha tree