From 4a16b4b3e997a74a124464360525b0aa5b3f46a1 Mon Sep 17 00:00:00 2001 From: pbastide Date: Mon, 17 Apr 2017 15:33:54 +0200 Subject: [PATCH 1/6] pPCA before PhyloEM --- .../multivariate_estimation_SUN_rBM.R | 18 +++++-- simulations_study/multivariate_exploitation.R | 52 ++++++++++++++----- 2 files changed, 54 insertions(+), 16 deletions(-) diff --git a/simulations_study/multivariate_estimation_SUN_rBM.R b/simulations_study/multivariate_estimation_SUN_rBM.R index 10cd95d..f432c48 100644 --- a/simulations_study/multivariate_estimation_SUN_rBM.R +++ b/simulations_study/multivariate_estimation_SUN_rBM.R @@ -65,7 +65,17 @@ nbrSim <- length(simlist) ###################### ## Estimation Function ###################### -estimations_several_K <- function(X){ +estimations_several_K <- function(X, pPCA = FALSE){ + if (pPCA){ + Y_data <- t(X$Y_data) + if (anyNA(Y_data)) return(list(sim = X, res = NULL)) + rownames(Y_data) <- trees[[paste0(X$ntaxa)]]$tip.label + ## Do a pPCA + Y_data_pPCA <- phytools::phyl.pca(trees[[paste0(X$ntaxa)]], Y_data) + Y_data <- t(Y_data_pPCA$S) + } else { + Y_data <- X$Y_data + } alpha_grid <- find_grid_alpha(trees[[paste0(X$ntaxa)]], nbr_alpha = 10, factor_up_alpha = 2, @@ -193,7 +203,7 @@ registerDoParallel(cl) time_alpha_gird_fav <- system.time( simestimations_fav <- foreach(i = simlist[favorables], .packages = reqpckg) %dopar% { - estimations_several_K(i) + estimations_several_K(i, pPCA=TRUE) } ) # Stop the cluster (parallel) @@ -213,9 +223,9 @@ registerDoParallel(cl) ## Parallelized estimations time_alpha_gird_unfav <- system.time( - simestimations_unfav <- foreach(i = simlist[!favorables][1:3], .packages = reqpckg) %dopar% + simestimations_unfav <- foreach(i = simlist[!favorables], .packages = reqpckg) %dopar% { - estimations_several_K(i) + estimations_several_K(i, pPCA=TRUE) } ) # Stop the cluster (parallel) diff --git a/simulations_study/multivariate_exploitation.R b/simulations_study/multivariate_exploitation.R index 9c5bf8e..309871a 100644 --- a/simulations_study/multivariate_exploitation.R +++ b/simulations_study/multivariate_exploitation.R @@ -1172,41 +1172,69 @@ datestamp_day <- "2016-11-28" ak <- "_both_alpha" # ak <- "" load(paste0(saveresultfile, "_aranged", ak, "-", datestamp_day, ".RData")) +results_summary_K_select_true$pPCA <- FALSE +summary_scores_K_select_true$pPCA <- FALSE ## l1OU saveresultfile = "../Results/Simulations_Multivariate/multivariate_estimations_l1ou" datestamp_day <- "2016-10-31" load(paste0(saveresultfile, "_aranged", "-", datestamp_day, ".RData")) +results_summary_l1ou$pPCA <- FALSE +summary_scores_l1ou$pPCA <- FALSE ## Merging library(plyr) results_summary_l1ou$alpha_known <- FALSE -results_summary_K_select_true <- rbind.fill(results_summary_K_select_true, - results_summary_l1ou) -rm(results_summary_l1ou) +results_summary_K_select_true_all <- rbind.fill(results_summary_K_select_true, + results_summary_l1ou) +rm(results_summary_l1ou, results_summary_K_select_true) summary_scores_l1ou$alpha_known <- FALSE -summary_scores_K_select_true <- rbind.fill(summary_scores_K_select_true, - summary_scores_l1ou) -rm(summary_scores_l1ou) +summary_scores_K_select_true_all <- rbind.fill(summary_scores_K_select_true, + summary_scores_l1ou) +rm(summary_scores_l1ou, summary_scores_K_select_true) ## l1OU - PCA saveresultfile = "../Results/Simulations_Multivariate/multivariate_estimations_l1ou" datestamp_day <- "2017-01-24" load(paste0(saveresultfile, "_aranged", "-", datestamp_day, ".RData")) -levels(results_summary_l1ou$K_type) <- "l1ou_PCA" -levels(summary_scores_l1ou$K_type) <- "l1ou_PCA" +# levels(results_summary_l1ou$K_type) <- "l1ou_PCA" +# levels(summary_scores_l1ou$K_type) <- "l1ou_PCA" +results_summary_l1ou$pPCA <- TRUE +summary_scores_l1ou$pPCA <- TRUE results_summary_l1ou$alpha_known <- FALSE -results_summary_K_select_true <- rbind.fill(results_summary_K_select_true, - results_summary_l1ou) +results_summary_K_select_true_all <- rbind.fill(results_summary_K_select_true_all, + results_summary_l1ou) rm(results_summary_l1ou) summary_scores_l1ou$alpha_known <- FALSE -summary_scores_K_select_true <- rbind.fill(summary_scores_K_select_true, - summary_scores_l1ou) +summary_scores_K_select_true_all <- rbind.fill(summary_scores_K_select_true_all, + summary_scores_l1ou) rm(summary_scores_l1ou) +## EM - PCA +saveresultfile = "../Results/Simulations_Multivariate/multivariate_estimations_SUN_rBM_pPCA" +datestamp_day <- "2017-04-18" +load(paste0(saveresultfile, "_aranged", "-", datestamp_day, ".RData")) +results_summary_K_select_true$pPCA <- TRUE +summary_scores_K_select_true$pPCA <- TRUE +results_summary_K_select_true$alpha_known <- FALSE +summary_scores_K_select_true$alpha_known <- FALSE + +results_summary_K_select_true_all <- rbind.fill(results_summary_K_select_true_all, + results_summary_K_select_true) +rm(results_summary_K_select_true) + +summary_scores_K_select_true_all <- rbind.fill(summary_scores_K_select_true_all, + summary_scores_K_select_true) +rm(summary_scores_K_select_true) + +## Correct names +results_summary_K_select_true <- results_summary_K_select_true_all +summary_scores_K_select_true <- summary_scores_K_select_true_all +rm(results_summary_K_select_true_all, summary_scores_K_select_true_all) + saveresultfile = "../Results/Simulations_Multivariate/multivariate_estimations_SUN_rBM_l1ou_PCA" datestamp_day <- format(Sys.time(), "%Y-%m-%d") save.image(paste0(saveresultfile, "_aranged", "-", datestamp_day, ".RData")) From f239eb0a55c638c705e1b5cdd31f15a2a7e8b81b Mon Sep 17 00:00:00 2001 From: Paul Bastide Date: Wed, 19 Apr 2017 09:50:36 +0200 Subject: [PATCH 2/6] Typo pPCA. --- simulations_study/multivariate_estimation_SUN_rBM.R | 2 +- simulations_study/multivariate_exploitation.R | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/simulations_study/multivariate_estimation_SUN_rBM.R b/simulations_study/multivariate_estimation_SUN_rBM.R index f432c48..71b036e 100644 --- a/simulations_study/multivariate_estimation_SUN_rBM.R +++ b/simulations_study/multivariate_estimation_SUN_rBM.R @@ -84,7 +84,7 @@ estimations_several_K <- function(X, pPCA = FALSE){ log_transform = TRUE) time_SUN <- system.time( res <- PhyloEM(phylo = trees[[paste0(X$ntaxa)]], - Y_data = X$Y_data, + Y_data = Y_data, process = "scOU", K_max = max(K_try[[paste0(X$ntaxa)]]) + 5, random.root = TRUE, diff --git a/simulations_study/multivariate_exploitation.R b/simulations_study/multivariate_exploitation.R index 309871a..545760e 100644 --- a/simulations_study/multivariate_exploitation.R +++ b/simulations_study/multivariate_exploitation.R @@ -1215,7 +1215,7 @@ rm(summary_scores_l1ou) ## EM - PCA saveresultfile = "../Results/Simulations_Multivariate/multivariate_estimations_SUN_rBM_pPCA" -datestamp_day <- "2017-04-18" +datestamp_day <- "2017-04-20" load(paste0(saveresultfile, "_aranged", "-", datestamp_day, ".RData")) results_summary_K_select_true$pPCA <- TRUE summary_scores_K_select_true$pPCA <- TRUE From 7fd45a9dd519b26a0f9db65d580b567e7024d3ab Mon Sep 17 00:00:00 2001 From: Paul Bastide Date: Tue, 25 Apr 2017 19:54:25 +0200 Subject: [PATCH 3/6] Add init.c for registration of native functions. --- NEWS.md | 2 ++ R/E_step.R | 2 +- src/init.c | 28 ++++++++++++++++++++++++++++ 3 files changed, 31 insertions(+), 1 deletion(-) create mode 100644 src/init.c diff --git a/NEWS.md b/NEWS.md index b9e166a..0284319 100644 --- a/NEWS.md +++ b/NEWS.md @@ -5,6 +5,8 @@ * Plotting missing values correctly in plot.PhyloEM. * Bug fixes in plotting PhyloEM object when p = 1. * When p=1 and nbr_alpha is not missing, do not switch to estimated mode for alpha. +* Technical: + * registration of c++ code to comply with R 3.4 new standards. # PhylogeneticEM 1.0.0 Initial Release. diff --git a/R/E_step.R b/R/E_step.R index cfafe47..02ef980 100644 --- a/R/E_step.R +++ b/R/E_step.R @@ -966,7 +966,7 @@ compute_log_likelihood.simple.nomissing.BM <- function(phylo, Y_data, sim, ## Upward_downward ############################################################################### -#' @useDynLib PhylogeneticEM +#' @useDynLib PhylogeneticEM, .registration=TRUE #' @importFrom Rcpp evalCpp # @import RcppArmadillo diff --git a/src/init.c b/src/init.c new file mode 100644 index 0000000..3cd6226 --- /dev/null +++ b/src/init.c @@ -0,0 +1,28 @@ +#include +#include +#include // for NULL +#include + +/* FIXME: + Check these declarations against the C/Fortran source code. + */ + +/* .Call calls */ +extern SEXP PhylogeneticEM_log_likelihood_BM(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP); +extern SEXP PhylogeneticEM_log_likelihood_OU(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP); +extern SEXP PhylogeneticEM_upward_downward_BM(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP); +extern SEXP PhylogeneticEM_upward_downward_OU(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP); + +static const R_CallMethodDef CallEntries[] = { + {"PhylogeneticEM_log_likelihood_BM", (DL_FUNC) &PhylogeneticEM_log_likelihood_BM, 6}, + {"PhylogeneticEM_log_likelihood_OU", (DL_FUNC) &PhylogeneticEM_log_likelihood_OU, 7}, + {"PhylogeneticEM_upward_downward_BM", (DL_FUNC) &PhylogeneticEM_upward_downward_BM, 6}, + {"PhylogeneticEM_upward_downward_OU", (DL_FUNC) &PhylogeneticEM_upward_downward_OU, 7}, + {NULL, NULL, 0} +}; + +void R_init_PhylogeneticEM(DllInfo *dll) +{ + R_registerRoutines(dll, NULL, CallEntries, NULL, NULL); + R_useDynamicSymbols(dll, FALSE); +} \ No newline at end of file From e67b2de1be06b8d4cfd52d55ea16ab61a504d92a Mon Sep 17 00:00:00 2001 From: Paul Bastide Date: Fri, 28 Apr 2017 16:16:11 +0200 Subject: [PATCH 4/6] Added arguments `label_font` and `axis_las` to plot. --- NAMESPACE | 2 +- NEWS.md | 3 +++ R/plot_functions.R | 15 +++++++++++++-- man/plot.PhyloEM.Rd | 9 +++++++-- 4 files changed, 24 insertions(+), 5 deletions(-) diff --git a/NAMESPACE b/NAMESPACE index 74735bd..d602af8 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -105,4 +105,4 @@ importFrom(utils,capture.output) importFrom(utils,combn) importFrom(utils,setTxtProgressBar) importFrom(utils,txtProgressBar) -useDynLib(PhylogeneticEM) +useDynLib(PhylogeneticEM, .registration=TRUE) diff --git a/NEWS.md b/NEWS.md index 0284319..36b45dc 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,6 +1,9 @@ # PhylogeneticEM 1.0.0.900 * Minor Changes: * Impose a maximum value for alpha in find_grid_alpha to respect machine max double. +* New Features + * added argument `label_font` to `plot` function to control the label font. + * added argument `axis_las` to `plot` function to control the axis las. * Bug fixes: * Plotting missing values correctly in plot.PhyloEM. * Bug fixes in plotting PhyloEM object when p = 1. diff --git a/R/plot_functions.R b/R/plot_functions.R index 770129e..e06583f 100644 --- a/R/plot_functions.R +++ b/R/plot_functions.R @@ -240,9 +240,11 @@ edgelabels_home <- function (text, edge, adj = c(0.5, 0.5), frame = "rect", #' @param show.tip.label whether to show the tip labels. Default to FALSE. #' @param label_cex if \code{show.tip.label=TRUE}, the size of the labels. Default #' to 0.5. +#' @param label_font if \code{show.tip.label=TRUE}, the font of the labels (see \link{par}). #' @param label_offset if \code{show.tip.label=TRUE}, the size of the offset between #' the tree and the labels. Default to 0. #' @param axis_cex cex for the label values of the plot. Default to 0.7. +#' @param axis_las las for the label values of the plot. Default to 0 (see \link{par}). #' @param edge.width width of the edge. Default to 1. #' @param margin_plot vector giving the margin to around the plot. #' Default to \code{c(0, 0, 0, 0)}. @@ -283,8 +285,10 @@ plot.PhyloEM <- function(x, alpha_border = 70, show.tip.label = FALSE, label_cex = 0.5, + label_font = 1, label_offset = 0, axis_cex = 0.7, + axis_las = 0, edge.width = 1, margin_plot = NULL, gray_scale = FALSE, @@ -357,6 +361,7 @@ plot.PhyloEM <- function(x, value_in_box = value_in_box, shifts_cex = shifts_cex, axis_cex = axis_cex, + axis_las = axis_las, margin_plot = margin_plot, color_shifts_regimes = color_shifts_regimes, # shifts_regimes = shifts_regimes, @@ -366,6 +371,7 @@ plot.PhyloEM <- function(x, ancestral_cex = ancestral_cex, ancestral_pch = ancestral_pch, label_cex = label_cex, + label_font = label_font, show.tip.label = show.tip.label, # underscore = underscore, # label.offset = label.offset, @@ -392,6 +398,7 @@ plot.data.process.actual <- function(Y.state, phylo, params, value_in_box = TRUE, shifts_cex = 1, axis_cex = 0.7, + axis_las = 0, margin_plot = NULL, color_shifts_regimes = FALSE, # shifts_regimes = NULL, @@ -401,6 +408,7 @@ plot.data.process.actual <- function(Y.state, phylo, params, ancestral_cex = 2, ancestral_pch = 19, label_cex = 1, + label_font = 1, show.tip.label = FALSE, underscore = FALSE, label.offset = 0, @@ -556,7 +564,9 @@ plot.data.process.actual <- function(Y.state, phylo, params, axis(1, at = pos_last_tip + eccart_g + range(Y.plot, na.rm = TRUE), labels = round(range(Y.state[t, ], na.rm = TRUE), digits = 2), pos = y.lim.min + ntaxa/15, - cex.axis = axis_cex, padj = -0.5) + cex.axis = axis_cex, + # padj = -0.5, + las = axis_las) # segments(pos_last_tip + eccart_g, y.lim.min + ntaxa/15, # pos_last_tip + eccart_g + unit, y.lim.min + ntaxa/15, # lwd = 2) @@ -588,7 +598,8 @@ plot.data.process.actual <- function(Y.state, phylo, params, if (!exists("color_characters_regimes")) color_characters_regimes <- color_characters text(x.lim.max.data, lastPP$yy[1:ntaxa], phylo$tip.label, cex = label_cex, pos = 4, - col = as.vector(color_characters_regimes)) + col = as.vector(color_characters_regimes), + font = label_font) } } ## Ancestral states diff --git a/man/plot.PhyloEM.Rd b/man/plot.PhyloEM.Rd index c6fd937..987ed59 100644 --- a/man/plot.PhyloEM.Rd +++ b/man/plot.PhyloEM.Rd @@ -12,8 +12,9 @@ shifts_cex = 0.6, shifts_bg = "chocolate4", root_bg = "chocolate4", shifts_adj = 0, root_adj = 1, color_shifts_regimes = FALSE, regime_boxes = FALSE, alpha_border = 70, show.tip.label = FALSE, - label_cex = 0.5, label_offset = 0, axis_cex = 0.7, edge.width = 1, - margin_plot = NULL, gray_scale = FALSE, root.edge = TRUE, ...) + label_cex = 0.5, label_font = 1, label_offset = 0, axis_cex = 0.7, + axis_las = 0, edge.width = 1, margin_plot = NULL, gray_scale = FALSE, + root.edge = TRUE, ...) } \arguments{ \item{x}{an object of class \code{PhyloEM}, result of function @@ -92,11 +93,15 @@ the border of the box. Default to 70.} \item{label_cex}{if \code{show.tip.label=TRUE}, the size of the labels. Default to 0.5.} +\item{label_font}{if \code{show.tip.label=TRUE}, the font of the labels (see \link{par}).} + \item{label_offset}{if \code{show.tip.label=TRUE}, the size of the offset between the tree and the labels. Default to 0.} \item{axis_cex}{cex for the label values of the plot. Default to 0.7.} +\item{axis_las}{las for the label values of the plot. Default to 0 (see \link{par}).} + \item{edge.width}{width of the edge. Default to 1.} \item{margin_plot}{vector giving the margin to around the plot. From 01055185d6b02f46ac58f51167dc9cddb634b57f Mon Sep 17 00:00:00 2001 From: Paul Bastide Date: Mon, 1 May 2017 10:31:19 +0200 Subject: [PATCH 5/6] Getting ready for CRAN 1.0.1 --- DESCRIPTION | 2 +- NEWS.md | 2 +- cran-comments.md | 4 +--- 3 files changed, 3 insertions(+), 5 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index a863b8c..3f18257 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: PhylogeneticEM Title: Automatic Shift Detection using a Phylogenetic EM -Version: 1.0.0.9000 +Version: 1.0.1 Authors@R: c( person("Paul", "Bastide", email = "paul.bastide@m4x.org", role = c("aut", "cre")), person("Mahendra", "Mariadassou", role = "ctb")) diff --git a/NEWS.md b/NEWS.md index 36b45dc..cdc80a8 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,4 +1,4 @@ -# PhylogeneticEM 1.0.0.900 +# PhylogeneticEM 1.0.1 * Minor Changes: * Impose a maximum value for alpha in find_grid_alpha to respect machine max double. * New Features diff --git a/cran-comments.md b/cran-comments.md index 95d8bcb..29c1fd0 100644 --- a/cran-comments.md +++ b/cran-comments.md @@ -1,5 +1,5 @@ ## Test environments -* ubuntu 12.04 (on travis-ci), R 3.3.2 +* ubuntu 12.04 and 16.04 (on travis-ci and local), R 3.4.0 * win-builder (devel and release) ## R CMD check results @@ -10,8 +10,6 @@ There was one NOTE: * checking CRAN incoming feasibility ... NOTE Maintainer: 'Paul Bastide ' -New submission - Possibly mis-spelled words in DESCRIPTION: OU (9:48) Ornstein–Uhlenbeck (9:29) From 683c8faea8b1370de452c02c06c92876f7e34486 Mon Sep 17 00:00:00 2001 From: Paul Bastide Date: Mon, 1 May 2017 11:27:14 +0200 Subject: [PATCH 6/6] remove NEWS.md from .Rbuildignore --- .Rbuildignore | 1 - 1 file changed, 1 deletion(-) diff --git a/.Rbuildignore b/.Rbuildignore index 5937f22..a9169cd 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -15,5 +15,4 @@ ^docs$ ^data-raw$ ^README\.md$ -^NEWS\.md$ ^cran-comments\.md$