add final convergence checking of diffsel

lib/diffsel.ml

@@ -31,6 +31,32 @@ let diffsel ~(phy_n:nucleotide_phylip workflow) ~(tree: _ workflow) ~(w_every:in
     )
   ]
 
+let check_conv run_diffsel : text_file directory workflow =
+  let env = docker_image ~account:"carinerey" ~name:"r_basics" ~tag:"07162018" () in
+  let script = tmp // "DiffselMCMCConvergenceAnalysis.Rmd" in
+  let trace = run_diffsel / selector["myrun.trace"] in
+  let out = dest // "out.html" in
+  let nb_new_iterations = dest // "new_iterations.txt.txt" in
+  workflow ~descr:"convergence_detection.DiffselMCMCConvergenceAnalysis" [
+    docker env (
+      and_list [
+        mkdir_p tmp ;
+        mkdir_p dest ;
+        cd tmp ;
+        cmd "cp" [ file_dump (string Scripts.diffselMCMCConvergenceAnalysis) ; script] ;
+        cmd "Rscript" [
+          string "-e" ;
+          string {|"rmarkdown::render(\"DiffselMCMCConvergenceAnalysis.Rmd\",|} ;
+          string {|params=list(set_trace1=\"|} ;
+          dep trace ;
+          string {|\"))"|};
+        ] ;
+        cmd "cp" [string "DiffselMCMCConvergenceAnalysis.html" ; ident out] ;
+        cmd "cp" [string "new_iterations.txt" ; ident nb_new_iterations]
+        ]
+    )
+  ]
+
 let selector run_diffsel : text_file workflow =
   let env = docker_image ~account:"vlanore" ~name:"diffsel" ~tag:"v1.0" () in
   let package = tmp // "diffsel_script_utils.py" in
@@ -46,7 +72,6 @@ let selector run_diffsel : text_file workflow =
       and_list [
         mkdir_p tmp ;
         cd tmp ;
-        (* cmd "ls" [dep run_diffsel]; (\* FIXME: PV ça sert à quoi ? *\) *)
 
         cmd "cp" [dep_ali; tmp_ali]; (* required dep to link the file in the env *)
         cmd "cp" [dep_tree; tmp_tree]; (* required dep to link the file in the env *)

lib/diffsel.mli

@@ -13,3 +13,7 @@ val diffsel :
 val selector :
   [`diffsel] directory workflow ->
   text_file workflow
+
+val check_conv :
+  [`diffsel] directory workflow ->
+  text_file directory workflow

lib/pipeline.ml

@@ -140,8 +140,7 @@ let derive_sim ~tree_dir ~trees ~profile_fn ~preview ~use_concat ~ns ~no_Ne ~no_
 let repo_of_detection_result res =
   let det_meth_prefix = Convergence_detection.meth_string_of_result res in
   Repo.[
-    (
-      match res with
+    [ match res with
       | `Pcoc w -> item ["pcoc.results.tsv"] (Pcoc.results w)
       | `Pcoc_gamma w -> item ["pcoc_gamma.results.tsv"] (Pcoc.results w)
       | `Pcoc_C60 w -> item ["pcoc_C60.results.tsv"] (Pcoc.results w)
@@ -152,8 +151,8 @@ let repo_of_detection_result res =
       | `Topological_WAG w -> item ["Topological_WAG.results.tsv"] (Topological.results w)
       | `Tdg09 w -> item ["Tdg09.results.tsv"] (Tamuri.results w)
       | `Multinomial w -> item ["Multinomial.results.tsv"] (Multinomial.results w)
-    ) ;
-    (
+    ] ;
+    [
       match res with
       | `Pcoc w -> item ["raw_results"]  w
       | `Pcoc_gamma w -> item ["raw_results"] w
@@ -165,8 +164,11 @@ let repo_of_detection_result res =
       | `Topological_WAG w -> item ["raw_results"] w
       | `Tdg09 w -> item ["raw_results"] w
       | `Multinomial w -> item ["raw_results"] w
-    ) ;
-  ]
+    ] ;
+     match res with
+       | `Diffsel w -> [item ["chain_convergence_checking.html"] ((Diffsel.check_conv w) / selector ["out.html"])]
+       | _ -> []
+  ] |> List.concat
   |> Repo.shift det_meth_prefix
   |> Repo.shift "Detection_tools"
 

lib/scripts/diffselMCMCConvergenceAnalysis.Rmd

+---
+title: "Diffsel MCMC convergence study"
+author: "Reviewphiltrans (Bastien Boussau)"
+output: html_document
+params:
+  set_trace1: "myrun.trace"
+---
+
+```{r setup, include=FALSE}
+knitr::opts_chunk$set(echo = TRUE)
+
+```
+
+
+# Study of the convergence of diffsel, using 1 chain.
+We use the files .trace output by diffsel.
+
+## Reading data
+
+```{r, echo = F}
+trace1 = gsub(" ", "",params$set_trace1)
+```
+
+```{r, echo = T}
+library(coda)
+d1<-read.table(trace1, h=T, comment.char="&")[,1:7]
+```
+## Summary
+
+```{r, echo = T}
+md1 <- mcmc(d1)
+print(summary(md1))
+```
+
+## Autocorrelation plots
+
+```{r}
+a<-autocorr.plot(d1, auto.layout = F)
+
+```
+
+
+## Effective sizes
+
+### ESS
+
+```{r}
+d1ess <- apply(d1, 2, effectiveSize)
+plot(d1ess, main= "ESS of chain 1")
+```
+
+## Trace plots
+
+```{r}
+for (i in 1:length(d1[1,])) {
+  plot(d1[,1], t="l", ylab=colnames(d1)[i])
+}
+
+```
+
+## Convergence in quantiles (Raftery and Lewis) 
+
+We require that estimation of the 0.025 quantile should be precise at the 1% level.
+
+```{r}
+md1 <- mcmc(d1)
+print(summary(md1))
+test = raftery.diag(md1, r=0.01)
+
+print(test)
+
+required_it = as.numeric(test$resmatrix[2])
+run_it = dim(md1)[1]
+
+new_iterations = required_it - run_it
+new_iterations_str = if (new_iterations > 0) {as.character(new_iterations)} else {"0"}
+writeLines(new_iterations_str, "new_iterations.txt")
+```
+
+`r run_it ` iterations were already run.
+
+`r required_it ` iterations are needed.
+
+`r new_iterations` additionnal iterations must be run.
+
+