full multinomial benchmark

lib/multinomial.ml

@@ -14,9 +14,17 @@ let multinomial ~(tree_sc:_ pworkflow) ~(faa:aminoacid_fasta pworkflow) : text_f
     ]
   ]
 
-let%pworkflow multinomial_asymptotic_lrt ~(tree_sc:_ pworkflow) ~(faa:aminoacid_fasta pworkflow) (* : text_file pworkflow *) =
+let%pworkflow multinomial_ocaml_implementation ~meth ~(tree_sc:_ pworkflow) ~(faa:aminoacid_fasta pworkflow) (* : text_file pworkflow *) =
   let open Phylogenetics in
   let open Convdet in
+  let module MT = Multinomial_test in
+  let meth = [%param meth] in
+  let test = match meth with
+    | `Asymptotic_LRT -> MT.LRT.asymptotic_test
+    | `Simulation_LRT -> MT.LRT.simulation_test ~sample_size:10_000
+    | `Asymptotic_sparse -> MT.Sparse.asymptotic_test
+    | `Simulation_sparse -> MT.Sparse.simulation_test ~sample_size:10_000
+  in
   let fold_leaves (root : _ Tree.t) ~init ~f =
     let open Tree in
     let rec node acc branch_data n =
@@ -56,12 +64,10 @@ let%pworkflow multinomial_asymptotic_lrt ~(tree_sc:_ pworkflow) ~(faa:aminoacid_
       )
   in
   let site i =
-    let module MT = Multinomial_test in
-
     let c0 = (counts seqs0 i :> int array) in
     let c1 = (counts seqs1 i :> int array) in
     let d = MT.data ~x1:c0 ~x2:c1 in
-    let r = MT.LRT.asymptotic_test d in
+    let r = test d in
     r._T_, r.pvalue
   in
   let n = Alignment.ncols alignment in
@@ -73,3 +79,8 @@ let%pworkflow multinomial_asymptotic_lrt ~(tree_sc:_ pworkflow) ~(faa:aminoacid_
   in
   let header = "Sites\t1MinusLRT\tLikelihoodRatio" in
   Out_channel.write_lines [%dest] (header :: res_lines)
+
+let multinomial_asymptotic_lrt ~tree_sc ~faa = multinomial_ocaml_implementation ~meth:`Asymptotic_LRT ~tree_sc ~faa
+let multinomial_simulation_lrt ~tree_sc ~faa = multinomial_ocaml_implementation ~meth:`Simulation_LRT ~tree_sc ~faa
+let multinomial_asymptotic_sparse ~tree_sc ~faa = multinomial_ocaml_implementation ~meth:`Asymptotic_sparse ~tree_sc ~faa
+let multinomial_simulation_sparse ~tree_sc ~faa = multinomial_ocaml_implementation ~meth:`Simulation_sparse ~tree_sc ~faa

lib/pipeline2.ml

@@ -154,6 +154,21 @@ let multinomial_asymptotic_lrt d =
     ~tree_sc:(tree d)
     ~faa:(amino_acid_alignment d)
 
+let multinomial_asymptotic_sparse d =
+  Multinomial.multinomial_asymptotic_sparse
+    ~tree_sc:(tree d)
+    ~faa:(amino_acid_alignment d)
+
+let multinomial_simulation_lrt d =
+  Multinomial.multinomial_simulation_lrt
+    ~tree_sc:(tree d)
+    ~faa:(amino_acid_alignment d)
+
+let multinomial_simulation_sparse d =
+  Multinomial.multinomial_simulation_sparse
+    ~tree_sc:(tree d)
+    ~faa:(amino_acid_alignment d)
+
 let diffseldsparse ?pi ?shiftprob ?eps d =
   Diffseldsparse.diffseldsparse
     ?pi ?shiftprob ?eps
@@ -188,9 +203,12 @@ let oracle d =
 let multinomial_benchmark d =
   Utils.rp_curve
     ~oracle:(oracle d)
-    ~labels:[|"LRT"|]
+    ~labels:[|"LRT";"LRTsim";"sparse";"sparse_sim"|]
     ~results:[
       multinomial_asymptotic_lrt d ;
+      multinomial_simulation_lrt d ;
+      multinomial_asymptotic_sparse d ;
+      multinomial_simulation_sparse d ;
     ]
 
 let benchmark d =