Commit 4b3d32db authored by Philippe Veber's avatar Philippe Veber
Browse files

tk/Mutsel_sim_cpg: compute only rate vectors instead of rate matrices

quadratic coefficient decreases from 1.671e-05 to 1.212e-05.

> df <- data.frame(n = c(500,1000,1300,2000), t = c(3.62,12.77,20.77,49.07)) ; fit <- lm(t ~ I(n ^ 2), data = df) ; summary(fit)

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

Residuals:
       1        2        3        4
 0.05496  0.11786 -0.24227  0.06946

Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.360e-01  1.594e-01   3.362   0.0782 .
I(n^2)      1.212e-05  7.145e-08 169.576 3.48e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.2005 on 2 degrees of freedom
Multiple R-squared:  0.9999,	Adjusted R-squared:  0.9999
F-statistic: 2.876e+04 on 1 and 2 DF,  p-value: 3.477e-05
parent 8f248d0d
...@@ -114,6 +114,33 @@ module Evolution_model = struct ...@@ -114,6 +114,33 @@ module Evolution_model = struct
| None -> 0. | None -> 0.
) )
let rate_vector { nucleotide_rates ; omega ; scaled_fitness = _F_ ; gBGC ; rate_CpG ; context_CpG; _ } p =
let nuc_rates = (nucleotide_rates :> Nucleotide.matrix) in
NSCodon.Table.init (fun q ->
if NSCodon.equal p q then 0.
else
match NSCodon.neighbours p q with
| Some (pos, x_a, x_b) ->
let rate_CpG = if rate_increase_CpG (pos, x_a, x_b) context_CpG (inspect_codon p) then rate_CpG else 1. in
let _B_ = match Nucleotide.(inspect x_a, inspect x_b) with
| (A | T), (C | G) -> gBGC
| (C | G), (A | T) -> -. gBGC
| _ -> 0.
in
let selection_coefficient =
_B_ +.
if NSCodon.synonym p q then 0.
else
let aa_p = NSCodon.aa_of_codon p in
let aa_q = NSCodon.aa_of_codon q in
_F_.Amino_acid.%(aa_q) -. _F_.Amino_acid.%(aa_p)
in
let p_fix = fixation_probability selection_coefficient in
let q_ab = nuc_rates.Nucleotide.%{x_a, x_b} in
q_ab *. omega *. p_fix *. rate_CpG
| None -> 0.
)
let stationary_distribution p = NSCodon_rate_matrix.stationary_distribution (rate_matrix p) let stationary_distribution p = NSCodon_rate_matrix.stationary_distribution (rate_matrix p)
let nucleotides_of_codons codons = let nucleotides_of_codons codons =
...@@ -177,15 +204,13 @@ module Make(BI: Simulator.Branch_info) = struct ...@@ -177,15 +204,13 @@ module Make(BI: Simulator.Branch_info) = struct
y y
let sequence_gillespie_direct tree ~root ~param = let sequence_gillespie_direct tree ~root ~param =
let codon_rates = memo Evolution_model.rate_matrix in let codon_rates = memo (fun (p, s) -> Evolution_model.rate_vector p s) in
Tree.propagate tree ~init:root ~node:Fn.const ~leaf:Fn.const ~branch:(fun seq b -> Tree.propagate tree ~init:root ~node:Fn.const ~leaf:Fn.const ~branch:(fun seq b ->
let rec loop state remaining_time = let rec loop state remaining_time =
let rate_matrices = Array.mapi state ~f:(fun i _ -> codon_rates (param state i b)) in let n = Array.length state in
let rates = let rates =
Array.mapi rate_matrices ~f:(fun i mat -> Array.init n ~f:(fun i ->
NSCodon.Table.init (fun m -> codon_rates (param state i b, state.(i))
if NSCodon.equal m state.(i) then 0. else (mat :> NSCodon.matrix).NSCodon.%{state.(i), m}
)
) in ) in
let pos_rates = Array.map rates ~f:(fun r -> Owl.Stats.sum (r :> float array)) in let pos_rates = Array.map rates ~f:(fun r -> Owl.Stats.sum (r :> float array)) in
let total_rate = Array.reduce_exn pos_rates ~f:( +. ) in let total_rate = Array.reduce_exn pos_rates ~f:( +. ) 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