Simulator: produce histogram of fitness values

lib/pipeline2.ml

@@ -128,6 +128,7 @@ let rec nucleotide_alignment = function
     let tree = tree sim in
     let fitness_profiles = Workflow.input profiles in
     Simulator.simulator ~branch_factor ~n_ha ~n_h0 ~ne_s ~gBGC ~tree ~seed ~fitness_profiles ()
+    |> fst
 
 let amino_acid_alignment d = Bppsuite.fna2faa (nucleotide_alignment d)
 

lib/simulator.ml

@@ -68,11 +68,32 @@ let%pworkflow simulator ?branch_factor ?seed ~n_h0 ~n_ha ~ne_s:(ne_s0, ne_s1) ~g
   let species_name =
     Phylogenetics.Tree.leaves tree
     |> List.map ~f:(fun { name } -> Option.value_exn name) in
-  Out_channel.with_file [%dest] ~f:(fun oc ->
-      List.iter2_exn species_name ali ~f:(fun description sequence ->
-          fprintf oc ">%s\n%s\n" description (sequence :> string)
+  let save_alignment fn =
+    Out_channel.with_file fn ~f:(fun oc ->
+        List.iter2_exn species_name ali ~f:(fun description sequence ->
+            fprintf oc ">%s\n%s\n" description (sequence :> string)
+          )
+      )
+  in
+  let save_fitness_histogram dest =
+    let data =
+      Array.fold (Array.append h0_params ha_params) ~init:[] ~f:(fun acc (p,q) ->
+          (p.scaled_fitness :> float array) :: (q.scaled_fitness :> float array) :: acc
         )
-    )
+      |> Array.concat
+    in
+    OCamlR_grDevices.pdf dest ;
+    ignore (OCamlR_graphics.hist data) ;
+    OCamlR_grDevices.dev_off ()
+  in
+  Core.Unix.mkdir_p [%dest] ;
+  save_alignment (Filename.concat [%dest] "simulation.fa") ;
+  save_fitness_histogram (Filename.concat [%dest] "fitness_histogram.pdf")
+
+let simulator ?branch_factor ?seed ~n_h0 ~n_ha ~ne_s ~gBGC ~tree ~fitness_profiles () =
+  let dir = simulator ?branch_factor ?seed ~n_h0 ~n_ha ~ne_s ~gBGC ~tree ~fitness_profiles () in
+  Bistro.Workflow.select dir ["simulation.fa"],
+  Bistro.Workflow.select dir ["fitness_histogram.pdf"]
 
 let%pworkflow pair_tree ~branch_length1 ~branch_length2 ~npairs =
   let open Phylogenetics in

lib/simulator.mli

@@ -11,7 +11,7 @@ val simulator :
   tree:nhx pworkflow ->
   fitness_profiles:#text_file pworkflow ->
   unit ->
-  nucleotide_fasta pworkflow
+  nucleotide_fasta pworkflow * pdf pworkflow
 
 val pair_tree :
   branch_length1:float ->