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

Simulator: fixed serious bug

the parameters used for root distribution and site evolution used to
be different, because of unwanted repeated sampling. That's what you
get using impure functions.
parent 7df238ca
......@@ -22,14 +22,12 @@ let%pworkflow simulator ?branch_factor ?seed ~n_h0 ~n_ha ~ne_s:(ne_s0, ne_s1) ~g
|> Convdet.Amino_acid.vector_of_array_exn
|> rescale_fitness beta
in
let random_h0_param () =
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)
let h0_params = Array.init n_h0 ~f:(fun _ ->
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
(p, q)
)
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
......@@ -37,21 +35,29 @@ let%pworkflow simulator ?branch_factor ?seed ~n_h0 ~n_ha ~ne_s:(ne_s0, ne_s1) ~g
| None -> assert false
| Some (_, i) -> i
in
let rec random_ha_param () =
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 ; gBGC = gBGC0 } in
let q = { base_param with scaled_fitness = q ; gBGC = gBGC1 } in
(function
| 0 -> p
| 1 -> q
| _ -> assert false)
let ha_params = Array.init n_ha ~f:(fun _ ->
let rec loop () =
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 loop ()
else
let p = { base_param with scaled_fitness = p ; gBGC = gBGC0 } in
let q = { base_param with scaled_fitness = q ; gBGC = gBGC1 } in
(p, q)
in
loop ()
)
in
let params i cond =
let p, q =
if i < n_h0 then h0_params.(i)
else ha_params.(i - n_h0)
in
match cond with
| 0 -> p
| 1 -> q
| _ -> assert false
in
let params i =
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 ->
Convdet.Mutsel.stationary_distribution (params i root_condition)
......
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