Commit ea76b505 authored by Philippe Veber's avatar Philippe Veber
Browse files

Simulator: update wrt convdet

parent 5ce44fd1
......@@ -7,10 +7,12 @@ let%pworkflow simulator ?branch_factor ?seed ~n_h0 ~n_ha ~ne_s ~tree ~fitness_pr
let ne_s = [%param ne_s] in
let branch_factor = [%param branch_factor] in
let tree = Convdet.Simulator.tree_from_file ?alpha:branch_factor [%path tree] in
let fitness_profiles = Convdet.Fitness_profiles.read [%path fitness_profiles] in
let fitness_profiles =
Convdet.Profile_tsv.(read [%path fitness_profiles] |> to_fitness)
|> Array.map ~f:(Array.map ~f:(( *. ) ne_s)) in
let base_param =
let p = Convdet.Mutsel.random_param ~alpha_nucleotide:10. ~alpha_fitness:0.1 in
{ p with ne_s ; omega = 1. }
{ p with omega = 1. }
in
let random_profile () =
Random.int (Array.length fitness_profiles)
......@@ -19,7 +21,7 @@ let%pworkflow simulator ?branch_factor ?seed ~n_h0 ~n_ha ~ne_s ~tree ~fitness_pr
in
let random_h0_param () =
let p = random_profile () in
{ base_param with fitness = p }
{ base_param with scaled_fitness = p }
in
let most_probable_aa (pref : float Convdet.Amino_acid.vector) =
let arr = Array.mapi (pref :> float array) ~f:(fun i x -> x, i) in
......@@ -32,8 +34,8 @@ let%pworkflow simulator ?branch_factor ?seed ~n_h0 ~n_ha ~ne_s ~tree ~fitness_pr
let q = random_profile () in
if p = q || most_probable_aa p = most_probable_aa q then random_ha_param ()
else
let p = { base_param with fitness = p } in
let q = { base_param with fitness = q } in
let p = { base_param with scaled_fitness = p } in
let q = { base_param with scaled_fitness = q } in
(function
| 0 -> p
| 1 -> q
......@@ -43,7 +45,12 @@ let%pworkflow simulator ?branch_factor ?seed ~n_h0 ~n_ha ~ne_s ~tree ~fitness_pr
if i < n_h0 then Fn.const (random_h0_param ())
else random_ha_param () in
let root_condition = Convdet.Simulator.Mutsel.root_condition tree in
let ali = Convdet.Simulator.Mutsel.alignment tree ~root_condition ~n:(n_h0 + n_ha) params in
let root_dists = Array.init (n_h0 + n_ha) ~f:(fun i ->
Convdet.Mutsel.stationary_distribution (params i root_condition)
)
in
let root = Convdet.Simulator.Mutsel.hmm0 ~len:(n_h0 + n_ha) ~dist:(Array.get root_dists) in
let ali = Convdet.Simulator.Mutsel.alignment tree ~root params in
let species_name =
Phylogenetics.Tree.leaves tree
|> List.map ~f:(fun { name } -> Option.value_exn name) in
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment