📣 An issue occured with the embedded container registry on October 25 2021, between 10:30 and 12:10 (UTC+2). Any persisting issues should be reported to CC-IN2P3 Support. 🐛

Commit 7fa09ac2 authored by Philippe Veber's avatar Philippe Veber
Browse files

whitespacecommit

parent 1be6cf2f
...@@ -6,7 +6,7 @@ let () = ...@@ -6,7 +6,7 @@ let () =
let stationary_distribution = Amino_acid.random_profile 0.5 in let stationary_distribution = Amino_acid.random_profile 0.5 in
let exchangeability_matrix = Rate_matrix.Amino_acid.make (fun aa_i aa_j-> let exchangeability_matrix = Rate_matrix.Amino_acid.make (fun aa_i aa_j->
let i, j = Amino_acid.(to_int aa_i, to_int aa_j) in let i, j = Amino_acid.(to_int aa_i, to_int aa_j) in
if i <= j then float_of_int i *. 0.01 +. float_of_int j *. 0.001 if i <= j then float_of_int i *. 0.01 +. float_of_int j *. 0.001
else float_of_int j *. 0.01 +. float_of_int i *. 0.001 else float_of_int j *. 0.01 +. float_of_int i *. 0.001
) in ) in
let p = { let p = {
...@@ -16,11 +16,11 @@ let () = ...@@ -16,11 +16,11 @@ let () =
} in } in
let t = 100. in let t = 100. in
let transition_matrix_diag = Evolution_model.transition_probability_matrix p in let transition_matrix_diag = Evolution_model.transition_probability_matrix p in
let transition_matrix_exp = let transition_matrix_exp =
let m = Evolution_model.rate_matrix p in let m = Evolution_model.rate_matrix p in
fun t -> Amino_acid.Matrix.(expm (scal_mul t m)) fun t -> Amino_acid.Matrix.(expm (scal_mul t m))
in in
Core.Command.run (Bench.make_command [ Core.Command.run (Bench.make_command [
Bench.Test.create ~name:"diag" (fun () -> transition_matrix_diag t) ; Bench.Test.create ~name:"diag" (fun () -> transition_matrix_diag t) ;
Bench.Test.create ~name:"exp" (fun () -> transition_matrix_exp t) ; Bench.Test.create ~name:"exp" (fun () -> transition_matrix_exp t) ;
]) ])
\ No newline at end of file
...@@ -22,8 +22,8 @@ module Evolution_model = struct ...@@ -22,8 +22,8 @@ module Evolution_model = struct
) )
let transition_probability_matrix p = let transition_probability_matrix p =
let module V = Amino_acid.Vector in let module V = Amino_acid.Vector in
let module M = Amino_acid.Matrix in let module M = Amino_acid.Matrix in
let m = rate_matrix p in let m = rate_matrix p in
let sqrt_pi = V.map p.stationary_distribution ~f:Float.sqrt in let sqrt_pi = V.map p.stationary_distribution ~f:Float.sqrt in
let diag_pi = M.diagm sqrt_pi in let diag_pi = M.diagm sqrt_pi in
...@@ -41,7 +41,7 @@ module Evolution_model = struct ...@@ -41,7 +41,7 @@ module Evolution_model = struct
let stationary_distribution = Amino_acid.random_profile 0.5 in let stationary_distribution = Amino_acid.random_profile 0.5 in
let exchangeability_matrix = Rate_matrix.Amino_acid.make (fun aa_i aa_j-> let exchangeability_matrix = Rate_matrix.Amino_acid.make (fun aa_i aa_j->
let i, j = Amino_acid.(to_int aa_i, to_int aa_j) in let i, j = Amino_acid.(to_int aa_i, to_int aa_j) in
if i <= j then float_of_int i *. 0.01 +. float_of_int j *. 0.001 if i <= j then float_of_int i *. 0.01 +. float_of_int j *. 0.001
else float_of_int j *. 0.01 +. float_of_int i *. 0.001 else float_of_int j *. 0.01 +. float_of_int i *. 0.001
) in ) in
let p = { let p = {
...@@ -50,12 +50,12 @@ module Evolution_model = struct ...@@ -50,12 +50,12 @@ module Evolution_model = struct
scale=1.; scale=1.;
} in } in
let t = 100. in let t = 100. in
let diag_exp_matrix = transition_probability_matrix p t in let diag_exp_matrix = transition_probability_matrix p t in
let m = rate_matrix p in let m = rate_matrix p in
let exp_matrix = Amino_acid.Matrix.(expm (scal_mul t m)) in let exp_matrix = Amino_acid.Matrix.(expm (scal_mul t m)) in
Amino_acid.Matrix.robust_equal ~tol:1e-10 diag_exp_matrix exp_matrix Amino_acid.Matrix.robust_equal ~tol:1e-10 diag_exp_matrix exp_matrix
let%test "Matrix exponential through diagonalisation matches naive implementation" = let%test "Matrix exponential through diagonalisation matches naive implementation" =
test_diagonal_matrix_exponential () test_diagonal_matrix_exponential ()
end 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