CamlPaml.PhyloLik.node_posterior: avoid unnecessary computation of al…

…l betas when computing posterior distribution for the root. This commit approximately doubles the speed of PhyloCSF, since it was calling node_posterior on the root to compute the ancestral composition score.

lib/CamlPaml/PhyloLik.ml

@@ -112,7 +112,10 @@ let ensure_beta x =
 				inter.{a} <- bp.{a} *. (Gsl.Blas.dot (Gsl.Matrix.row ss a) xas)
 
 			for b = 0 to k-1 do
-				for a = 0 to k-1 do ps_colb.{a} <- ps.{a,b}
+				for a = 0 to k-1 do
+					(* idea: instead of this loop it'd probably be a little faster
+					         to transpose ps *)
+					Bigarray.Array1.unsafe_set ps_colb a (Bigarray.Array2.unsafe_get ps a b)
 				x.beta.{i,b} <- Gsl.Blas.dot inter ps_colb
 		x.have_beta <- true
 		
@@ -123,17 +126,16 @@ let likelihood x =
 
 let node_posterior inferred i =
 	if inferred.workspace.generation <> inferred.my_generation then failwith "CamlPaml.PhyloLik.node_posterior: invalidated workspace"
-	ensure_beta inferred
+	ensure_alpha inferred
 	let k = snd (Gsl.Matrix.dims inferred.alpha)
-	if inferred.z = 0. then
+	if inferred.z = 0. then (* the data were impossible by the model *)
 		Array.make k 0.
+	else if i < T.leaves inferred.tree then (* leaf: usually certain *)
+		Gsl.Vector.to_array (alpha_get inferred i)
 	else
-		let nleaves = T.leaves inferred.tree
-
-		if i < nleaves then
-			Gsl.Vector.to_array (alpha_get inferred i)
-		else
-			Array.init k (fun x -> (alpha_get inferred i).{x} *. inferred.beta.{i,x} /. inferred.z)
+		(* calculate the betas (except for the root, whose betas get prefilled in prepare) *)
+		if i <> T.root inferred.tree then ensure_beta inferred
+		Array.init k (fun x -> (alpha_get inferred i).{x} *. inferred.beta.{i,x} /. inferred.z)
 
 let add_branch_posteriors ?(weight=1.0) inferred branch ecounts =
 	if inferred.workspace.generation <> inferred.my_generation then failwith "CamlPaml.PhyloLik.add_branch_posterior: invalidated workspace"

src/test.ml

@@ -25,16 +25,18 @@ let run_PhyloCSF species params =
     | _ -> raise Exit
 
 let talAA () =
-  let ans = run_PhyloCSF "12flies" "../PhyloCSF_Examples/tal-AA.fa"
+  let ans = run_PhyloCSF "12flies" "../PhyloCSF_Examples/tal-AA.fa --ancComp"
   print_endline (String.join "\t" ans)
   List.nth ans 1 $hould # equal "score(decibans)"
   float_of_string (List.nth ans 2) $hould # be # within (297.62,297.63)
+  float_of_string (List.nth ans 3) $hould # be # within (48.25,48.26)
 
 let aldh2_ex5_out () =
-  let ans = run_PhyloCSF "29mammals" "../PhyloCSF_Examples/ALDH2.exon5.fa"
+  let ans = run_PhyloCSF "29mammals" "../PhyloCSF_Examples/ALDH2.exon5.fa --ancComp"
   print_endline (String.join "\t" ans)
   List.nth ans 1 $hould # equal "score(decibans)"
   float_of_string (List.nth ans 2) $hould # be # within (-178.93,-178.92)
+  float_of_string (List.nth ans 3) $hould # be # within (-38.29,-38.28)
 
 let aldh2_ex5_in () =
   let abbreviate = (try ignore (Sys.getenv "SKIP_SLOW"); true with Not_found -> false)