Merge branch 'develop'

.Rbuildignore

@@ -15,5 +15,4 @@
 ^docs$
 ^data-raw$
 ^README\.md$
-^NEWS\.md$
 ^cran-comments\.md$

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"))

NAMESPACE

@@ -105,4 +105,4 @@ importFrom(utils,capture.output)
 importFrom(utils,combn)
 importFrom(utils,setTxtProgressBar)
 importFrom(utils,txtProgressBar)
-useDynLib(PhylogeneticEM)
+useDynLib(PhylogeneticEM, .registration=TRUE)

NEWS.md

@@ -1,10 +1,15 @@
-# 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
+  * 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.
   * 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.

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
 

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

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 <paul.bastide@m4x.org>'
 
-New submission
-
 Possibly mis-spelled words in DESCRIPTION:
   OU (9:48)
   Ornstein–Uhlenbeck (9:29)

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,
@@ -74,7 +84,7 @@ estimations_several_K <- function(X){
                                 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,
@@ -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)

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-20"
+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"))

src/init.c

@@ -0,0 +1,28 @@
+#include <R.h>
+#include <Rinternals.h>
+#include <stdlib.h> // for NULL
+#include <R_ext/Rdynload.h>
+
+/* 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);
+}