TDG09 : faster matrix exponentiation (hopefully ?)

lib/tk/tdg09.ml

@@ -22,8 +22,13 @@ module Evolution_model = struct
       )
   let transition_probability_matrix p =
     let m = rate_matrix p in
+    let (d_vec, transform_matrix) = Amino_acid.Matrix.diagonalize m in
+    let rev_trans_matrix = Amino_acid.Matrix.inverse transform_matrix in
     fun t ->
-      Amino_acid.Matrix.(expm (scal_mul t m))
+      let exp_matrix = Amino_acid.Vector.(exp (scal_mul t d_vec))
+                       |> Amino_acid.Matrix.diagm in
+      Amino_acid.Matrix.(dot (dot transform_matrix exp_matrix) rev_trans_matrix)
+      (* Amino_acid.Matrix.(expm (scal_mul t m)) *)
 end
 
 let choose_aa p =
@@ -533,9 +538,9 @@ module Implementation_check = struct
     ignore (
       let values = Array.map results ~f:(fun (_,_,_,lrt) -> lrt.pvalue) in
       OCamlR_graphics.hist
-         ~main:title
-         ~xlab:"p"
-         ~breaks:(`n 20) values :> OCamlR_graphics.hist) ;
+        ~main:title
+        ~xlab:"p"
+        ~breaks:(`n 20) values :> OCamlR_graphics.hist) ;
     OCamlR_grDevices.dev_off ()
 
   let render_stat_histogram ~title ~df results dest =
@@ -543,10 +548,10 @@ module Implementation_check = struct
     ignore (
       let values = Array.map results ~f:(fun (_,_,_,lrt) -> lrt._D_) in
       OCamlR_graphics.hist
-         ~main:title
-         ~xlab:"D"
-         ~freq:false
-         ~breaks:(`n 20) values :> OCamlR_graphics.hist) ;
+        ~main:title
+        ~xlab:"D"
+        ~freq:false
+        ~breaks:(`n 20) values :> OCamlR_graphics.hist) ;
     let x = Array.init 1_000 ~f:(fun i -> float i /. 10.) in
     let y = Array.map x ~f:(Gsl.Randist.chisq_pdf ~nu:df) in
     OCamlR_graphics.lines ~x ~y () ;