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

tk/Mutsel_simulator_cpg: use gen-like iterator to encode outdated positions

parent 7411d86f
......@@ -193,6 +193,11 @@ module Make(BI: Simulator.Branch_info) = struct
Gsl.Randist.(discrete rng (discrete_preproc v))
|> NSCodon.of_int_exn
let cpg_update ~n ~pos f =
f pos ;
if pos - 1 >= 0 then f (pos - 1) ;
if pos + 1 < n then f (pos + 1)
let sequence_gillespie_direct tree ~root ~param =
Tree.propagate tree ~init:root ~node:Fn.const ~leaf:Fn.const ~branch:(fun seq b ->
let state = Array.copy seq in
......@@ -208,13 +213,8 @@ module Make(BI: Simulator.Branch_info) = struct
else
let pos = Discrete_pd.draw pos_rates rng in
let next_letter = symbol_sample (rates.(pos) :> float array) in
let outdated_positions =
[pos]
|> (fun l -> if pos - 1 >= 0 then (pos - 1) :: l else l)
|> (fun l -> if pos + 1 < n then (pos + 1) :: l else l)
in
state.(pos) <- next_letter ;
List.iter outdated_positions ~f:(fun pos ->
cpg_update ~n ~pos (fun pos ->
rates.(pos) <- rate pos ;
Discrete_pd.update pos_rates pos (pos_rate pos)
) ;
......
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