Commit 29e841f0 authored by Philippe Veber's avatar Philippe Veber
Browse files

Simulator: allow variable NeS and gBGC

parent 2923e9ad
......@@ -32,8 +32,8 @@ type dataset =
profiles : string ;
n_h0 : int ;
n_ha : int ;
ne_s : float ;
gBGC : float ;
ne_s : float * float ;
gBGC : float * float ;
seed : int ;
}
......@@ -56,7 +56,7 @@ let bppseqgen_simulation ~hyp ~tree ~profiles ~nb_sites ~seed =
seed ;
}
let convdet_simulation ?(branch_factor = 1.) ?(ne_s = 1.) ?(gBGC = 0.) ?(seed = 0) ~tree ~profiles ~n_h0 ~n_ha () =
let convdet_simulation ?(branch_factor = 1.) ?(ne_s = 1., 1.) ?(gBGC = 0., 0.) ?(seed = 0) ~tree ~profiles ~n_h0 ~n_ha () =
Convdet_simulation {
tree ;
profiles ;
......
open Core_kernel
let%pworkflow simulator ?branch_factor ?seed ~n_h0 ~n_ha ~ne_s ~gBGC ~tree ~fitness_profiles () =
let%pworkflow simulator ?branch_factor ?seed ~n_h0 ~n_ha ~ne_s:(ne_s0, ne_s1) ~gBGC:(gBGC0, gBGC1) ~tree ~fitness_profiles () =
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_s = [%param ne_s] in
let gBGC = [%param gBGC] in
let ne_s0 = [%param ne_s0] in
let ne_s1 = [%param ne_s1] in
let gBGC0 = [%param gBGC0] in
let gBGC1 = [%param gBGC1] 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.Profile_tsv.(read [%path fitness_profiles] |> to_fitness)
|> Array.map ~f:(Array.map ~f:(( *. ) ne_s)) in
let fitness_profiles = Convdet.Profile_tsv.(read [%path fitness_profiles] |> to_fitness) in
let rescale_fitness beta = Convdet.Amino_acid.vector_map ~f:(( *. ) beta) in
let base_param =
let p = Convdet.Mutsel.random_param ~alpha_nucleotide:10. ~alpha_fitness:0.1 in
{ p with omega = 1. ; gBGC }
{ p with omega = 1. }
in
let random_profile () =
let random_profile beta =
Random.int (Array.length fitness_profiles)
|> Array.get fitness_profiles
|> Convdet.Amino_acid.vector_of_array_exn
|> rescale_fitness beta
in
let random_h0_param () =
let p = random_profile () in
{ base_param with scaled_fitness = p }
let profile = random_profile 1. in
let p = { base_param with gBGC = gBGC0 ; scaled_fitness = rescale_fitness ne_s0 profile } in
let q = { base_param with gBGC = gBGC1 ; scaled_fitness = rescale_fitness ne_s1 profile } in
(function
| 0 -> p
| 1 -> q
| _ -> assert false)
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
......@@ -31,19 +38,19 @@ let%pworkflow simulator ?branch_factor ?seed ~n_h0 ~n_ha ~ne_s ~gBGC ~tree ~fitn
| Some (_, i) -> i
in
let rec random_ha_param () =
let p = random_profile () in
let q = random_profile () in
let p = random_profile ne_s0 in
let q = random_profile ne_s1 in
if p = q || most_probable_aa p = most_probable_aa q then random_ha_param ()
else
let p = { base_param with scaled_fitness = p } in
let q = { base_param with scaled_fitness = q } in
let p = { base_param with scaled_fitness = p ; gBGC = gBGC0 } in
let q = { base_param with scaled_fitness = q ; gBGC = gBGC1 } in
(function
| 0 -> p
| 1 -> q
| _ -> assert false)
in
let params i =
if i < n_h0 then Fn.const (random_h0_param ())
if i < n_h0 then random_h0_param ()
else random_ha_param () in
let root_condition = Convdet.Simulator.Mutsel.root_condition tree in
let root_dists = Array.init (n_h0 + n_ha) ~f:(fun i ->
......
......@@ -6,8 +6,8 @@ val simulator :
?seed:int ->
n_h0:int ->
n_ha:int ->
ne_s:float ->
gBGC:float ->
ne_s:float * float ->
gBGC:float * float ->
tree:nhx pworkflow ->
fitness_profiles:#text_file pworkflow ->
unit ->
......
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