Commit 7411d86f authored by Philippe Veber's avatar Philippe Veber
Browse files

tk/Mutsel_simulator_cpg: simulator with linear time complexity

now we can generate a million sites within a minute

> df <- data.frame(n = c(10000,30000,100000,300000), t = c(1.08,2.25,6.27,18.03)) ; fit <- lm(t ~ n, data = df) ; summary(fit)

Call:
lm(formula = t ~ n, data = df)

Residuals:
       1        2        3        4
 0.01851  0.01931 -0.05290  0.01509

Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.769e-01  2.997e-02   15.91  0.00393 **
n           5.846e-05  1.886e-07  310.00 1.04e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.04325 on 2 degrees of freedom
Multiple R-squared:      1,	Adjusted R-squared:      1
F-statistic: 9.61e+04 on 1 and 2 DF,  p-value: 1.041e-05
parent 3ca117c0
......@@ -200,12 +200,13 @@ module Make(BI: Simulator.Branch_info) = struct
let rate i = Evolution_model.rate_vector (param state i b) state.(i) in
let rates = Array.init n ~f:rate in
let pos_rate i = Owl.Stats.sum (rates.(i) :> float array) in
let pos_rates = Array.init n ~f:pos_rate in
let rec loop total_rate remaining_time =
let pos_rates = Discrete_pd.init n ~f:pos_rate in
let rec loop remaining_time =
let total_rate = Discrete_pd.total_weight pos_rates in
let tau = Owl.Stats.exponential_rvs ~lambda:total_rate in
if Float.(tau > remaining_time) then state
else
let pos = Gsl.Randist.(discrete rng (discrete_preproc pos_rates)) in
let pos = Discrete_pd.draw pos_rates rng in
let next_letter = symbol_sample (rates.(pos) :> float array) in
let outdated_positions =
[pos]
......@@ -213,19 +214,13 @@ module Make(BI: Simulator.Branch_info) = struct
|> (fun l -> if pos + 1 < n then (pos + 1) :: l else l)
in
state.(pos) <- next_letter ;
let delta_total_rate =
List.map outdated_positions ~f:(fun pos ->
let old = pos_rates.(pos) in
rates.(pos) <- rate pos ;
pos_rates.(pos) <- pos_rate pos ;
pos_rates.(pos) -. old
)
|> List.fold ~init:0. ~f:( +. )
in
loop (total_rate +. delta_total_rate) Float.(remaining_time - tau)
List.iter outdated_positions ~f:(fun pos ->
rates.(pos) <- rate pos ;
Discrete_pd.update pos_rates pos (pos_rate pos)
) ;
loop Float.(remaining_time - tau)
in
let total_rate = Array.reduce_exn pos_rates ~f:( +. ) in
loop total_rate (BI.length b)
loop (BI.length b)
)
end
......
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