In-house LMM x Gemma : fixes not positive-definite matrix

app/lmm_benchmark.ml

@@ -43,7 +43,7 @@ let orthomam_echolocation = {
     Orthomam.tree_of_db
       ~branch_length_unit:`Amino_acid
       ~convergent_species:Orthomam.species_with_echolocation
-      (Codepitk.Orthomam_db.make "omm")
+      (Codepitk.Orthomam_db.make "/home/louis/Data/omm")
   ) ;
   rooted = false ;
   branch_scale = 1. ;
@@ -60,12 +60,12 @@ let meth ?(requires_rooted_tree = false) f label =
   { f ; label ; requires_rooted_tree }
 
 let methods = Pipeline.[
-  meth tdg09 "tdg09" ;
-  meth pcoc "pcoc" ;
+  (* meth tdg09 "tdg09" ;
+  meth pcoc "pcoc" ; *)
   (* meth pcoc_v2 ~col:3 "pcoc v2" ; *)
   meth (gemma ~lmm_test:`Score ~relatedness_mode:`Standardized) "gemma" ;
   meth inhouse_lmm "LMM" ;
-  meth (lmm_with_gemma_matrix ~relatedness_mode:`Standardized) "LMM_Gemma";
+  meth (lmm_with_gemma_matrix ~relatedness_mode:`Centered) "LMM_Gemma";
   meth multinomial_asymptotic_lrt "multinomial" ;
   meth topological "topological" ~requires_rooted_tree:true ;
 ]
@@ -89,6 +89,7 @@ let benchmark_rds ?(seed = 42) { tree = t ; rooted ; ne_s ; branch_scale ; _ } =
 let () =
   let open Bistro_utils.Repo in
   let datasets = [ besnard2009 ; rubisco ; oneline_rodent ; orthomam_echolocation ] in
+  (* let datasets = [ besnard2009 ] in *)
   let repo =
     List.map datasets ~f:(fun d ->
         item [d.label ^ ".rds"] (benchmark_rds d)
@@ -96,3 +97,22 @@ let () =
   in
   let loggers = [ Bistro_utils.Console_logger.create () ] in
   build_main ~loggers ~np:4 ~mem:(`GB 4) repo ~outdir:"res"
+
+(* module Top = Bistro_utils.Toplevel_eval.Make(struct
+  let np = 10
+  let mem = 8 
+end)() *)
+(* 
+let () =
+  let q =
+    Pipeline.query ~seed:42 ~tree:(NHX besnard2009.tree) ~branch_scale:1. ~ne_s:(1., 1.)
+      ~profiles:"example/aa_fitness/263SelectedProfiles.tsv"
+      ~n_h0:90 ~n_ha:10 ()
+  in
+  (* Gemma.calculate_relatedness_matrix
+    ~mode:`Standardized
+    ~genotype:(Gemma.genotype_of_fasta (Pipeline.amino_acid_alignment q)) 
+    ~phenotype:(Gemma.phenotype_of_tree (Pipeline.tree ~branch_length_unit:`Amino_acid q)) *)
+    Pipeline.lmm_with_gemma_matrix q ~relatedness_mode:`Standardized
+  |> Top.path
+  |> print_endline *)
\ No newline at end of file

lib/r_console.ml

+open Core_kernel
+open Bistro
+
+type assignment =
+  | Path : string * _ path workflow -> assignment
+
+module type Toplevel = sig
+  val eval : 'a workflow -> 'a
+end
+
+let path var_id w = Path (var_id, w)
+
+let rconsole (module T : Toplevel) code assignments =
+  let paths = 
+    List.map assignments ~f:(fun (Path (_, w)) -> Workflow.path w)
+    |> Workflow.list
+    |> T.eval
+  in
+  print_endline code ;
+  List.iter2_exn assignments paths ~f:(fun (Path (var_id, _)) path ->
+      Printf.printf "%s <- '%s'\n" var_id path
+    )

lib/r_console.mli

+open Bistro
+
+type assignment
+
+module type Toplevel = sig
+    val eval : 'a workflow -> 'a
+end
+
+val path :
+    string ->
+    _ path workflow ->
+    assignment 
+
+val rconsole :
+    (module Toplevel) ->
+    string ->
+    assignment list ->
+    unit    

lib/tk/inhouse_lmm.ml

@@ -109,8 +109,8 @@ let phenotypes_of_tree t =
    correlation matrix, provided in sparse representation, [n] being
    the dimension of the (square) matrix and [cors] the list of non-zero coefficients *)
 let transform_matrix n cors =
-  let sigma () = matrix_of_correlations n cors in
-  let mat_C = in_place (L.potrf ~up:false) (sigma ()) in
+  let sigma = matrix_of_correlations n cors in
+  let mat_C = in_place (L.potrf ~up:false) sigma in
   for i = 1 to L.Mat.dim1 mat_C do
     for j = i + 1 to L.Mat.dim2 mat_C do
       mat_C.{i, j} <- 0.
@@ -123,7 +123,11 @@ let load_relatedness_matrix filename =
     In_channel.read_lines filename
     |> List.mapi ~f:(fun i row -> 
         String.split row ~on:'\t'
-        |> List.mapi ~f:(fun j value -> i, j, (float_of_string value))
+        |> List.mapi ~f:(fun j value -> 
+          (* small addition on the diagonal to ensure eigen values are strictly positive *)
+          let m_ij = float_of_string value +. if i = j then 0.000001 else 0. in
+          i, j, m_ij
+      )
      ) in 
   transform_matrix (List.length rowwise_matrix) (List.concat rowwise_matrix)