CamlPaml: add test coverage of codepaths for non-reversible rate matr…

…ices

lib/CamlPaml/Makefile

@@ -5,7 +5,7 @@ lib:
 
 test:
 	ocamlbuild test.native
-	_build/test.native -verbose
+	_build/test.native
 
 
 install: lib

lib/CamlPaml/Q.ml

@@ -121,7 +121,7 @@ module Diag = struct
 		assert (n = n')
 		n
 
-	let of_Q ?(tol=1e-6) qm =
+	let of_Q ?(tol=1e-6) ?(force_complex=false) qm =
 		let qm = Gsl.Matrix.of_arrays qm
 		let l, s = Gsl.Eigen.nonsymmv ~protect:true (`M qm)
 		let s' = zinvm s
@@ -131,7 +131,7 @@ module Diag = struct
 		for i = 0 to n-1 do if not (check_real ~tol l.{i}) then rev := false
 
 		let eig =
-			if !rev then
+			if !rev && not force_complex then
 				`r { r_s = m_of_cm s; r_s' = m_of_cm s'; r_l = v_of_cv l }
 			else
 				`nr { nr_s = s; nr_s' = s'; nr_l = l }
@@ -144,13 +144,18 @@ module Diag = struct
 
 	let reversible = function
 		| { eig = `r _ } -> true
-		| _ -> false
+		| { eig = `nr { nr_l }; tol } ->
+			let rev = ref true
+			let n = Gsl.Vector_complex.length nr_l
+			for i = 0 to n-1 do if not (check_real ~tol:tol nr_l.{i}) then rev := false
+			!rev
 
 	let equilibrium q =
 		if not q.have_pi then
 			let eig = match q.eig with
 				| `r eig -> eig
-				| `nr _ -> failwith "CamlPaml.Q.equilibrium: non-reversible model"
+				| `nr { nr_s; nr_s'; nr_l } when reversible q -> { r_s = m_of_cm nr_s; r_s' = m_of_cm nr_s'; r_l = v_of_cv nr_l }
+				| _ -> failwith "CamlPaml.Q.equilibrium: non-reversible model"
 			let n = Gsl.Vector.length eig.r_l
 			let min_L = ref infinity
 			let min_Lp = ref (-1)

lib/CamlPaml/Q.mli

@@ -15,8 +15,11 @@ module Diag : sig
 	type t
 	
 	(** Diagonalize the rate matrix. The rate matrix needs not be reversible (the internal
-        representation uses complex arithmetic). *)
-	val of_Q : ?tol:float -> matrix -> t
+        representation can use complex arithmetic).
+
+    @param force_complex force internal use of complex arithmetic, even if the rate matrix
+                         is reversible. *)
+	val of_Q : ?tol:float -> ?force_complex:bool -> matrix -> t
 	val to_Q : t -> matrix
 	
 	(** return [n] for an [n-by-n] matrix *)

lib/CamlPaml/test.ml

@@ -72,17 +72,31 @@ module BeagleTests = struct
 		let qdiagtests = "Q.Diag" >::: [TestCase (fun () -> assert_bool "Q.Diag.reversible" (Q.Diag.reversible q));
 										TestCase (fun () -> assert_equal_vector ~msg:"Q.Diag.equilibrium" (Q.Diag.equilibrium q) [| 0.25; 0.25; 0.25; 0.25 |])]
 										
-		let likelihood =
-			TestCase
-				fun () ->
-					let m = PhyloModel.make t [| q |]
-					let ll = ref 0.
-					for i = 0 to String.length human - 1 do
-						let lf = Array.map (fun ch -> if Code.DNA.is ch then `Certain (Code.DNA.index ch) else `Marginalize) [| human.[i]; chimp.[i]; gorilla.[i] |]
-						ll := !ll +. log (PhyloLik.likelihood (PhyloModel.prepare_lik m lf))
-					assert_equal ~cmp:(cmp_float ~epsilon:0.001) ~printer:string_of_float ~msg:"mismatch" (-1574.63623) !ll
+		let likelihood which_q () =
+			let m = PhyloModel.make t [| which_q |]
+			let ll = ref 0.
+			for i = 0 to String.length human - 1 do
+				let lf = Array.map (fun ch -> if Code.DNA.is ch then `Certain (Code.DNA.index ch) else `Marginalize) [| human.[i]; chimp.[i]; gorilla.[i] |]
+				ll := !ll +. log (PhyloLik.likelihood (PhyloModel.prepare_lik m lf))
+			assert_equal ~cmp:(cmp_float ~epsilon:0.001) ~printer:string_of_float ~msg:"mismatch" (-1574.63623) !ll
+
+		let qc =
+			Q.Diag.scale
+				Q.Diag.of_Q ~force_complex:true [|
+						[| (-3.); 1.; 1.; 1.; |];
+						[| 1.; (-3.); 1.; 1.; |];
+						[| 1.; 1.; (-3.); 1.; |];
+						[| 1.; 1.; 1.; (-3.); |]
+					|]
+				(1. /. 3.)
+
+		let qcdiagtests = "Q.Diag (complex)" >::: [
+				TestCase (fun () -> assert_bool "Q.Diag.reversible" (Q.Diag.reversible qc));
+				TestCase (fun () -> assert_equal_vector ~msg:"Q.Diag.equilibrium" (Q.Diag.equilibrium qc) [| 0.25; 0.25; 0.25; 0.25 |])
+			]
 					
-		let tests = "TinyTest" >::: [ qdiagtests; "likelihood" >: likelihood ]
+		let tests = "TinyTest" >::: [ qdiagtests; "likelihood" >:: likelihood q;
+									  qcdiagtests; "likelihood (complex)" >:: likelihood qc ]
 
 	let tests = "BeagleTest" >::: [ TinyTest.tests ]