Mutsel CpG simulator: update make_param signature

Also adds some convenience functions for flat and random params

lib/tk/mutsel_simulator_cpg.ml

@@ -40,7 +40,7 @@ module Evolution_model = struct
     Amino_acid.Vector.init (fun _ -> 1. /. float Amino_acid.card)
     |> fitness_of_profile
 
-  let make_param state pos _branch =
+  let make_param ~gBGC ~scaled_fitness ~omega ~nucleotide_rates ~nucleotide_stat_dist ~rate_CpG state pos _branch =
     let context_CpG =
       if pos = 0 then
         context_CpG_of_codons state.(pos) ~rh:state.(pos+1) ()
@@ -49,19 +49,29 @@ module Evolution_model = struct
       else
         context_CpG_of_codons state.(pos) ~lh:state.(pos -1) ~rh:state.(pos + 1) ()
     in
-    let pi = Nucleotide.flat_profile () in
-    let rho = Linear_algebra.Lacaml.Vector.init 6 ~f:(fun _ -> 1. /. 6.) in
-    let nucleotide_rates = Rate_matrix.Nucleotide.gtr ~equilibrium_frequencies:pi ~transition_rates:rho in
     {
       nucleotide_rates;
-      nucleotide_stat_dist=pi;
-      omega=1.;
-      scaled_fitness = flat_fitness ();
-      gBGC = 0.;
-      rate_CpG = 0.;
+      nucleotide_stat_dist;
+      omega;
+      scaled_fitness;
+      gBGC;
+      rate_CpG;
       context_CpG;
     }
 
+  let flat_param () = 
+    let pi = Nucleotide.flat_profile () in
+    let rho = Linear_algebra.Lacaml.Vector.init 6 ~f:(fun _ -> 1. /. 6.) in
+    let nucleotide_rates = Rate_matrix.Nucleotide.gtr ~equilibrium_frequencies:pi ~transition_rates:rho in 
+    make_param ~gBGC:0. ~scaled_fitness:(flat_fitness ()) ~omega:1. ~nucleotide_rates ~nucleotide_stat_dist:pi ~rate_CpG:1. 
+
+  let random_param ~alpha_nuc ~alpha_fitness = 
+    let pi = Nucleotide.random_profile alpha_nuc in 
+    let rho = Utils.random_profile 6 in 
+    let nucleotide_rates = Rate_matrix.Nucleotide.gtr ~equilibrium_frequencies:pi ~transition_rates:rho in 
+    let scaled_fitness = Amino_acid.random_profile alpha_fitness |> fitness_of_profile in
+    make_param ~gBGC:0. ~scaled_fitness ~omega:1. ~nucleotide_rates ~nucleotide_stat_dist:pi ~rate_CpG:1.
+
   let rate_increase_CpG (pos, x_a, x_b) context =
     match Nucleotide.(pos, inspect x_a, inspect x_b), context with
     | (0, C, T), { codon = (_,G, _) ; _ }

lib/tk/mutsel_simulator_cpg.mli

@@ -23,7 +23,18 @@ module Evolution_model : sig
 
   val stationary_distribution:  param -> NSCodon.vector
 
-  val make_param: NSCodon.t array -> int -> 'a -> param
+  val make_param: 
+    gBGC:float -> 
+    scaled_fitness:Amino_acid.vector -> 
+    omega:float -> 
+    nucleotide_rates:Nucleotide.matrix ->  
+    nucleotide_stat_dist:Nucleotide.vector -> 
+    rate_CpG:float -> 
+    NSCodon.t array -> int -> 'a -> param
+
+  val flat_param: unit -> NSCodon.t array -> int -> 'a -> param
+
+  val random_param: alpha_nuc:float -> alpha_fitness:float -> NSCodon.t array -> int -> 'a -> param
 end
 
 module Make(BI: Simulator.Branch_info) : sig