Mutsel simulator stationary distribution

lib/tk/mutsel_simulator_cpg.ml

@@ -2,6 +2,7 @@ open Core_kernel
 open Phylogenetics
 
 module NSCodon = Codon.Universal_genetic_code.NS
+module NSCodon_rate_matrix = Rate_matrix.Make(NSCodon)
 
 module Evolution_model = struct
   type context_CpG = {
@@ -103,19 +104,7 @@ module Evolution_model = struct
         | None -> 0.
       )
 
-  let stationary_distribution p =
-    let pi = p.nucleotide_stat_dist in
-    NSCodon.Vector.init (fun codon ->
-        let n1, n2, n3 = NSCodon.nucleotides codon in
-        let aa = NSCodon.aa_of_codon codon in
-        let b n = match Nucleotide.inspect n with
-          | A | T -> -. p.gBGC
-          | C | G -> +. p.gBGC
-        in
-        Nucleotide.(pi.%(n1) *. pi.%(n2) *. pi.%(n3))
-        *. exp (p.scaled_fitness.Amino_acid.%(aa) +. (b n1 +. b n2 +. b n3) /. 2.)
-      )
-    |> NSCodon.Vector.normalize
-end
+  let stationary_distribution p = NSCodon_rate_matrix.stationary_distribution (rate_matrix p)
+end 
 
 module Make(BI: Simulator.Branch_info) = Simulator.Make(NSCodon)(Evolution_model)(BI)