From 88e0fc8d5ed54107cdb2411d23de718b2b6e6f0c Mon Sep 17 00:00:00 2001 From: Paul Bastide Date: Wed, 31 Aug 2022 11:07:06 +0200 Subject: [PATCH 01/10] Dev version --- DESCRIPTION | 2 +- NEWS.md | 2 ++ 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index f34ca1b..468d516 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: PhylogeneticEM Title: Automatic Shift Detection using a Phylogenetic EM -Version: 1.6.0 +Version: 1.6.0.9000 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 0b3b9bb..b7f5cce 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,5 @@ +# PhylogeneticEM 1.6.0.9000 + # PhylogeneticEM 1.6.0 * Deprecation fix: * Update code to comply with Matrix 1.4-2 new coding standards. From 122b071744ef58e17878847dc9807d3023f3b0eb Mon Sep 17 00:00:00 2001 From: Paul Bastide Date: Thu, 5 Oct 2023 09:22:57 +0200 Subject: [PATCH 02/10] Fix documentation issues. --- DESCRIPTION | 2 +- R/E_step.R | 24 +++++++++---------- R/M_step.R | 12 +++++----- R/estimateEM.R | 10 ++++---- R/generic_functions.R | 2 +- R/init_EM.R | 2 +- R/model_selection.R | 2 -- R/parsimonyNumber.R | 12 +++++----- R/partitionsNumber.R | 2 ++ R/shifts_manipulations.R | 3 +++ R/shutoff.R | 2 +- R/simulate.R | 18 +++++++------- man/allocate_subset_node.simulate.Rd | 6 ++--- man/check_data.Rd | 2 -- man/compute_E.simple.Rd | 8 ------- man/compute_expectations.BM.Rd | 4 ++-- man/compute_expectations.scOU.Rd | 4 ++-- man/compute_gauss_lasso.Rd | 4 ++-- man/compute_log_likelihood.simple.Rd | 6 ----- man/compute_mahalanobis_distance.simple.Rd | 6 ++--- man/compute_residuals.simple.Rd | 6 ++--- man/enlight.Rd | 4 ++-- man/equivalent_shifts_values.Rd | 4 ---- man/estimateEM.Rd | 14 +++++------ man/extract.variance_covariance.Rd | 4 ++-- man/find_shift_values.Rd | 2 -- man/go_back_to_original_process.Rd | 8 +++---- man/imputed_traits.Rd | 4 ++-- man/init.parsimonyCost.Rd | 4 ++-- man/init.simulate.OU.Rd | 8 +++---- man/init.simulate.StateAndExp.Rd | 8 +++---- man/is.in.ranges.Rd | 4 ++-- man/log_likelihood.Rd | 6 ++--- man/merge_params_independent.Rd | 2 +- man/model_selection.Rd | 4 ++-- man/penalty_pBIC.Rd | 3 --- man/{qr.solve_exact.Rd => qr_solve_exact.Rd} | 6 ++--- man/segmentation.OU.specialCase.lasso.Rd | 5 ---- ...segmentation.OU.specialCase.same_shifts.Rd | 5 ---- man/simul_process.Rd | 6 ++--- man/simulate_internal.Rd | 4 ++-- tests/testthat/test-utilities-rotations.R | 5 ++-- 42 files changed, 107 insertions(+), 140 deletions(-) rename man/{qr.solve_exact.Rd => qr_solve_exact.Rd} (89%) diff --git a/DESCRIPTION b/DESCRIPTION index 468d516..bb31afe 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -43,7 +43,7 @@ LinkingTo: License: GPL (>= 2) | file LICENSE Encoding: UTF-8 LazyData: true -RoxygenNote: 7.1.2 +RoxygenNote: 7.2.3 URL: https://github.com/pbastide/PhylogeneticEM BugReports: https://github.com/pbastide/PhylogeneticEM/issues VignetteBuilder: knitr diff --git a/R/E_step.R b/R/E_step.R index 95da816..419da49 100644 --- a/R/E_step.R +++ b/R/E_step.R @@ -33,10 +33,10 @@ #' \code{extract_simulate_internal} to extract the needed quantities from these objects. #' #' @param phylo Input tree. -#' @param Y_data : vector indicating the data at the tips -#' @param sim (list) : result of function \code{simulate} -#' @param Sigma : variance-covariance matrix, result of function \code{compute_variance_covariance} -#' @param Sigma_YY_inv : invert of the variance-covariance matrix of the data +# @param Y_data vector indicating the data at the tips +# @param sim (list) result of function \code{simulate} +# @param Sigma variance-covariance matrix, result of function \code{compute_variance_covariance} +# @param Sigma_YY_inv invert of the variance-covariance matrix of the data #' #' @return conditional_law_X (list) : list of conditional statistics : #' "expectation" : matrix of size p x (ntaxa+Nnode), with ntaxa @@ -314,11 +314,11 @@ compute_fixed_moments <- function(times_shared, ntaxa){ #' #' @param struct structural matrix of size (ntaxa+Nnode)*p, result #' of function \code{compute_variance_covariance} -#' @param what: sub-matrix to be extracted: +#' @param what sub-matrix to be extracted: #' "YY" : sub-matrix of tips (p*ntaxa first lines and columns) #' "YZ" : sub matrix tips x nodes (p*Nnode last rows and p*ntaxa first columns) #' "ZZ" : sub matrix of nodes (p*Nnode last rows and columns) -#' @param miss; missing values of Y_data +#' @param miss missing values of Y_data #' #' @return sub-matrix of variance covariance. #' @@ -902,9 +902,9 @@ compute_mean_variance.simple.nomissing.BM <- function (phylo, #' to extract the needed quantities from these objects. #' #' @param phylo Input tree. -#' @param Y_data : vector indicating the data at the tips. +#' @param Y_data_vec : vector indicating the data at the tips. #' @param sim (list) : result of function \code{simulate}. -#' @param Sigma_YY_inv : invert of the variance-covariance matrix of the data. +#' @param Sigma_YY_chol_inv : invert of the Cholesky variance-covariance matrix of the data. #' #' @keywords internal #' @@ -933,9 +933,9 @@ compute_residuals.simple <- function(phylo, Y_data_vec, sim, #' to extract the needed quantities from these objects. #' #' @param phylo Input tree. -#' @param Y_data : vector indicating the data at the tips. +#' @param Y_data_vec : vector indicating the data at the tips. #' @param sim (list) : result of function \code{simulate}. -#' @param Sigma_YY_inv : invert of the variance-covariance matrix of the data. +#' @param Sigma_YY_chol_inv : invert of the cholesky variance-covariance matrix of the data. #' #' @return squared Mahalanobis distance between data and mean at the tips. #' @@ -979,10 +979,8 @@ compute_mahalanobis_distance.simple.nomissing.BM <- function(phylo, Y_data, sim, #' \code{extract_simulate_internal} to extract the needed quantities from these objects. #' #' @param phylo Input tree. -#' @param Y_data : vector indicating the data at the tips -#' @param sim (list) : result of function \code{simulate} #' @param Sigma : variance-covariance matrix, result of function \code{compute_variance_covariance} -#' @param Sigma_YY_inv : invert of the variance-covariance matrix of the data +# @param Sigma_YY_inv : invert of the variance-covariance matrix of the data #' #' @return log likelihood of the data #' diff --git a/R/M_step.R b/R/M_step.R index 813a825..5d95d78 100644 --- a/R/M_step.R +++ b/R/M_step.R @@ -954,9 +954,9 @@ segmentation.BM <- function(nbr_of_shifts, costs0, diff_exp){ #' #' @param phylo a phylogenetic tree #' @param nbr_of_shifts Number of shifts on the phylogeny allowed -#' @param conditional_law_X moments of the conditional law of X given Y, result -#' of function \code{compute_M.OU.specialCase} -#' @param selection.strength the selection strength +# @param conditional_law_X moments of the conditional law of X given Y, result +# of function \code{compute_M.OU.specialCase} +# @param selection.strength the selection strength #' #' @return List containing : beta_0 : the optimal value at the root #' shifts : list containing the computed tau and delta @@ -1143,9 +1143,9 @@ compute_regression_matrices <- function(phylo, conditional_law_X, selection.stre #' This is the best move if keeping the previous shifts positions. #' #' @param phylo a phylogenetic tree -#' @param conditional_law_X moments of the conditional law of X given Y, result -#' of function \code{compute_M.OU.specialCase} -#' @param selection.strength the selection strength +# @param conditional_law_X moments of the conditional law of X given Y, result +# of function \code{compute_M.OU.specialCase} +# @param selection.strength the selection strength #' @param shifts_old the previous list of shifts (only position is used) #' #' @return List containing : beta_0 : the optimal value at the root diff --git a/R/estimateEM.R b/R/estimateEM.R index 6ca3586..163588f 100644 --- a/R/estimateEM.R +++ b/R/estimateEM.R @@ -2806,7 +2806,7 @@ estimateEM_wrapper_scratch <- function(phylo, Y_data, #' @param phylo a phylogenetic tree, class \code{\link[ape]{phylo}}. #' @param Y_data matrix of data at the tips (pxntaxa) #' @param check.tips.names (bool) whether to check the tips names or not -#' @param trait_correlation_threshold threshold for trait correlation. Default to 0.9. +# @param trait_correlation_threshold threshold for trait correlation. Default to 0.9. #' #' @return Y_data a re-ordered matrix of data (if necessary) #' @@ -2963,11 +2963,11 @@ return_to_original_order <- function(X, phy_o, phy_r){ #' on a rescaled tree, and gives back the equivalent parameters of the OU on #' the original process. #' -#' @param phy_original: the original phylogenetic tree -#' @param known.selection.strength: the known selection strength of the original +#' @param phy_original the original phylogenetic tree +#' @param known.selection.strength the known selection strength of the original #' OU. -#' @param sBM_variance: boolean. Is the root random ? -#' @param params: the inferred parameters of the BM on the re-scaled tree. +#' @param sBM_variance boolean. Is the root random ? +#' @param params the inferred parameters of the BM on the re-scaled tree. #' #' #' @return params_scOU the equivalent parameters of the OU on the original tree. diff --git a/R/generic_functions.R b/R/generic_functions.R index 110b896..8104782 100644 --- a/R/generic_functions.R +++ b/R/generic_functions.R @@ -860,7 +860,7 @@ split_params_independent <- function(params){ #' objects into one param object of dimension p #' The reverse operation is done by \code{split_params_independent} #' -#' @param params_split: a list of parameters +#' @param params_split a list of parameters #' #' @return A parameter object #' diff --git a/R/init_EM.R b/R/init_EM.R index 15208b2..f68ed67 100644 --- a/R/init_EM.R +++ b/R/init_EM.R @@ -1427,7 +1427,7 @@ find_independent_regression_vectors.gglasso <- function(Xkro, K, fit, root, p, g #' multiplication with a Cholesky decomposition of the variance) is included, rather #' than the intercept. #' -#' @param Yp (transformed) data +#' @param Ypt (transformed) data #' @param Xp (transformed) matrix of regression #' @param delta regression coefficients obtained with a lasso regression #' @param root the position of the root (intercept) in delta diff --git a/R/model_selection.R b/R/model_selection.R index 4120ac2..19fc4af 100644 --- a/R/model_selection.R +++ b/R/model_selection.R @@ -294,8 +294,6 @@ model_selection_BGH_leastsquares_raw <- function(res, ntaxa, C.LINselect, ...){ #' @param K the dimension of the model. #' @param model_complexity the complexity of the set of models with dimension K. #' @param ntaxa the number of tips. -#' @param C a constant, C > 1. Default is C = 1.1 -#' (as suggested in Baraud Giraud Huet (2009)) #' #' @return value of the penalty. #' diff --git a/R/parsimonyNumber.R b/R/parsimonyNumber.R index 61624e1..9654981 100644 --- a/R/parsimonyNumber.R +++ b/R/parsimonyNumber.R @@ -152,7 +152,7 @@ extract.parsimonyCost <- function(x, #' At a tip i in state k, the line-vector is initialized as follow : #' (1 - Ind(k=p)_{1<=p<=nclus})*Inf (where Inf * 0 = 0) #' -#' @param phylo a phylogenetic tree, class \code{\link[ape]{phylo}}. +#' @param phy a phylogenetic tree, class \code{\link[ape]{phylo}}. #' @param clusters the vector of the clusters of the tips. #' #' @return A (ntaxa + Nnode)x(nclus) matrix, with ntaxa first lines initialized as @@ -1401,8 +1401,8 @@ equivalent_shifts_edges <- function(phylo, #' the tips. Careful, only works for ULTRAMETRIC trees. #' #' @param phylo a phylogenetic tree, class \code{\link[ape]{phylo}}. -#' @param shifts a list of positions and values of original shifts. -#' @param beta_0 value of the original optimal value at the root. +# @param shifts a list of positions and values of original shifts. +# @param beta_0 value of the original optimal value at the root. #' @param eq_shifts_edges matrix (optional) result of function #' \code{\link{equivalent_shifts_edges}}. #' @param T_tree_ac matrix of incidence of the tree, result of function @@ -1447,7 +1447,7 @@ equivalent_shifts_values <- function(phylo, #' @param shifts_edges a vector of positions of shifts on the tree. #' @param T_tree_ac matrix of incidence of the tree, result of function #' \code{incidence.matrix}. -#' @param m_Y the vector of values of the means at the tips. +# @param m_Y the vector of values of the means at the tips. #' #' @return vector, with first entry the values at the root, and other entries the #' values of the shifts. @@ -1457,7 +1457,7 @@ equivalent_shifts_values <- function(phylo, find_shift_values <- function(shifts_edges, T_tree_ac, vect_Y, p, ntaxa){ pos_shifts <- as.vector(sapply(shifts_edges, function(z) p * (z-1) + 1:p)) mat <- cbind(t(matrix(rep(diag(1, p, p), ntaxa), nrow = p)), T_tree_ac[, pos_shifts]) # stationary case assumption used here - coefs <- try(qr.solve_exact(mat, vect_Y), silent = TRUE) + coefs <- try(qr_solve_exact(mat, vect_Y), silent = TRUE) if (inherits(coefs, "try-error")){ warning("Had a problem solving exactly the linear system.") return(rep(NA, length(shifts_edges) + 1)) @@ -1482,7 +1482,7 @@ find_shift_values <- function(shifts_edges, T_tree_ac, vect_Y, p, ntaxa){ #' #' @keywords internal ## -qr.solve_exact <- function (a, b, tol = 1e-07) { +qr_solve_exact <- function (a, b, tol = 1e-07) { if (!is.qr(a)) a <- qr(a, tol = tol) nc <- ncol(a$qr) diff --git a/R/partitionsNumber.R b/R/partitionsNumber.R index 0029bfe..e60cd1f 100644 --- a/R/partitionsNumber.R +++ b/R/partitionsNumber.R @@ -305,6 +305,8 @@ update.partitionsNumber.gen <- function(daughtersParams, ...){ prod.index <- function(X,Id){ if (0 %in% Id) return(0) # Indice hors limites return(prod(diag(X[,Id]))) + # if (any(Id == 0)) return(0) # Indice hors limites + # return(prod(X[1:length(Id) + (Id - 1) * nrow(X)])) } ## diff --git a/R/shifts_manipulations.R b/R/shifts_manipulations.R index aec2478..f866d82 100644 --- a/R/shifts_manipulations.R +++ b/R/shifts_manipulations.R @@ -487,6 +487,9 @@ update.compute_betas_from_shifts <- function(edgeNbr, ancestral, shifts, ...){ #' @return Vector of size (ntaxa + Nnode) of the optimal values at the tips #' of the tree. #' +#' @example +#' +#' #' @export #' #06/10/14 - Initial release diff --git a/R/shutoff.R b/R/shutoff.R index 9d61d97..0050fe9 100644 --- a/R/shutoff.R +++ b/R/shutoff.R @@ -213,7 +213,7 @@ shutoff.EM.OU.stationary.root_AND_shifts_at_nodes.half_life <- function(params_o #' @details #' This function is used to test the convergence of the algorithm. #' -#' @param params list of parameters with the correct structure +#' @param p list of parameters with the correct structure #' @param min list of minimum values for the parameters #' @param max list of maximum values for the parameters #' diff --git a/R/simulate.R b/R/simulate.R index 032a817..47de8fa 100644 --- a/R/simulate.R +++ b/R/simulate.R @@ -550,9 +550,9 @@ extract_simulate_internal <- function(paramSimu, #' \code{allocate_subset_node.simulate} slice the data correctly and allocate the right part. #' To be used in function UpdateDown. #' -#' @param node: node on which to slice -#' @param array: structure to be sliced -#' @param value: value to be attributed to the slice +#' @param node node on which to slice +#' @param array structure to be sliced +#' @param value value to be attributed to the slice #' #' @return array: array p x Nnode x 2 (BM), with slice corresponding to node filled with value #' @@ -575,8 +575,8 @@ subset_node.simulate <- function(node, array){ #' #' @description Function used in \code{\link{simulate}} for BM/OU initializations. #' -#' @param phy: Input tree. -#' @param p: dimension of the trait simulated +#' @param phy Input tree. +#' @param p dimension of the trait simulated #' @param root.state (list): state of the root, with: #' random : random state (TRUE) or deterministic state (FALSE) #' value.root : if deterministic, value of the character at the root @@ -642,8 +642,8 @@ init.simulate.BM <- function(phy, p, root.state, simulate_random, ...){ #' #' @description Function used in \code{\link{simulate}} for OU initialization. #' -#' @param phy: Input tree. -#' @param p: dimension of the trait simulated +#' @param phy Input tree. +#' @param p dimension of the trait simulated #' @param root.state (list): state of the root, with: #' random : random state (TRUE) or deterministic state (FALSE) #' value.root : if deterministic, value of the character at the root @@ -853,7 +853,7 @@ update.simulate.OUBM <- function(edgeNbr, ancestral, #' \code{compute_expectations.BM} use the matrix formulation to compute the #' expected values at all the nodes. #' -#' @param phylo: Input tree. +#' @param phylo Input tree. #' @param root.state (list): state of the root, with: #' random : random state (TRUE) or deterministic state (FALSE) #' value.root : if deterministic, value of the character at the root @@ -893,7 +893,7 @@ compute_expectations.BM <- function(phylo, root.state, shifts, U_tree = NULL){ #' \code{compute_expectations.scOU} use the matrix formulation to compute the #' expected values at all the nodes. Assumes a stationary root. #' -#' @param phylo: Input tree. +#' @param phylo Input tree. #' @param root.state (list): state of the root, with: #' random : random state (TRUE) or deterministic state (FALSE) #' value.root : if deterministic, value of the character at the root diff --git a/man/allocate_subset_node.simulate.Rd b/man/allocate_subset_node.simulate.Rd index 3d3daa4..e089c03 100644 --- a/man/allocate_subset_node.simulate.Rd +++ b/man/allocate_subset_node.simulate.Rd @@ -7,11 +7,11 @@ allocate_subset_node.simulate(node, array, value) } \arguments{ -\item{node:}{node on which to slice} +\item{node}{node on which to slice} -\item{array:}{structure to be sliced} +\item{array}{structure to be sliced} -\item{value:}{value to be attributed to the slice} +\item{value}{value to be attributed to the slice} } \value{ array: array p x Nnode x 2 (BM), with slice corresponding to node filled with value diff --git a/man/check_data.Rd b/man/check_data.Rd index 9211b94..a2ff576 100644 --- a/man/check_data.Rd +++ b/man/check_data.Rd @@ -12,8 +12,6 @@ check_data(phylo, Y_data, check.tips.names) \item{Y_data}{matrix of data at the tips (pxntaxa)} \item{check.tips.names}{(bool) whether to check the tips names or not} - -\item{trait_correlation_threshold}{threshold for trait correlation. Default to 0.9.} } \value{ Y_data a re-ordered matrix of data (if necessary) diff --git a/man/compute_E.simple.Rd b/man/compute_E.simple.Rd index d976ef4..ede7b62 100644 --- a/man/compute_E.simple.Rd +++ b/man/compute_E.simple.Rd @@ -22,14 +22,6 @@ compute_E.simple( } \arguments{ \item{phylo}{Input tree.} - -\item{Y_data}{: vector indicating the data at the tips} - -\item{sim}{(list) : result of function \code{simulate}} - -\item{Sigma}{: variance-covariance matrix, result of function \code{compute_variance_covariance}} - -\item{Sigma_YY_inv}{: invert of the variance-covariance matrix of the data} } \value{ conditional_law_X (list) : list of conditional statistics : diff --git a/man/compute_expectations.BM.Rd b/man/compute_expectations.BM.Rd index 2252c00..ed8fb46 100644 --- a/man/compute_expectations.BM.Rd +++ b/man/compute_expectations.BM.Rd @@ -7,6 +7,8 @@ compute_expectations.BM(phylo, root.state, shifts, U_tree = NULL) } \arguments{ +\item{phylo}{Input tree.} + \item{root.state}{(list): state of the root, with: random : random state (TRUE) or deterministic state (FALSE) value.root : if deterministic, value of the character at the root @@ -18,8 +20,6 @@ edges : vector of the K id of edges where the shifts are values : matrix p x K of values of the shifts on the edges (one column = one shift) relativeTimes : vector of dimension K of relative time of the shift from the parent node of edges} - -\item{phylo:}{Input tree.} } \value{ paramSimu: array p x Nnode x 2 (BM). For each trait t, 1 <= t <= p, diff --git a/man/compute_expectations.scOU.Rd b/man/compute_expectations.scOU.Rd index 24a5fc0..892d50c 100644 --- a/man/compute_expectations.scOU.Rd +++ b/man/compute_expectations.scOU.Rd @@ -14,6 +14,8 @@ compute_expectations.scOU( ) } \arguments{ +\item{phylo}{Input tree.} + \item{root.state}{(list): state of the root, with: random : random state (TRUE) or deterministic state (FALSE) value.root : if deterministic, value of the character at the root @@ -25,8 +27,6 @@ edges : vector of the K id of edges where the shifts are values : matrix p x K of values of the shifts on the edges (one column = one shift) relativeTimes : vector of dimension K of relative time of the shift from the parent node of edges} - -\item{phylo:}{Input tree.} } \value{ paramSimu: array p x Nnode x 2 (BM). For each trait t, 1 <= t <= p, diff --git a/man/compute_gauss_lasso.Rd b/man/compute_gauss_lasso.Rd index 3ffc31e..149604e 100644 --- a/man/compute_gauss_lasso.Rd +++ b/man/compute_gauss_lasso.Rd @@ -13,13 +13,13 @@ compute_gauss_lasso( ) } \arguments{ +\item{Ypt}{(transformed) data} + \item{Xp}{(transformed) matrix of regression} \item{delta}{regression coefficients obtained with a lasso regression} \item{root}{the position of the root (intercept) in delta} - -\item{Yp}{(transformed) data} } \value{ Named list, with "E0.gauss" the intercept (value at the root); diff --git a/man/compute_log_likelihood.simple.Rd b/man/compute_log_likelihood.simple.Rd index 5695cfc..c336555 100644 --- a/man/compute_log_likelihood.simple.Rd +++ b/man/compute_log_likelihood.simple.Rd @@ -19,13 +19,7 @@ compute_log_likelihood.simple( \arguments{ \item{phylo}{Input tree.} -\item{sim}{(list) : result of function \code{simulate}} - \item{Sigma}{: variance-covariance matrix, result of function \code{compute_variance_covariance}} - -\item{Y_data}{: vector indicating the data at the tips} - -\item{Sigma_YY_inv}{: invert of the variance-covariance matrix of the data} } \value{ log likelihood of the data diff --git a/man/compute_mahalanobis_distance.simple.Rd b/man/compute_mahalanobis_distance.simple.Rd index 6a2e852..971cd42 100644 --- a/man/compute_mahalanobis_distance.simple.Rd +++ b/man/compute_mahalanobis_distance.simple.Rd @@ -16,11 +16,11 @@ compute_mahalanobis_distance.simple( \arguments{ \item{phylo}{Input tree.} -\item{sim}{(list) : result of function \code{simulate}.} +\item{Y_data_vec}{: vector indicating the data at the tips.} -\item{Y_data}{: vector indicating the data at the tips.} +\item{sim}{(list) : result of function \code{simulate}.} -\item{Sigma_YY_inv}{: invert of the variance-covariance matrix of the data.} +\item{Sigma_YY_chol_inv}{: invert of the cholesky variance-covariance matrix of the data.} } \value{ squared Mahalanobis distance between data and mean at the tips. diff --git a/man/compute_residuals.simple.Rd b/man/compute_residuals.simple.Rd index 61543ec..1523c64 100644 --- a/man/compute_residuals.simple.Rd +++ b/man/compute_residuals.simple.Rd @@ -9,11 +9,11 @@ compute_residuals.simple(phylo, Y_data_vec, sim, Sigma_YY_chol_inv, miss) \arguments{ \item{phylo}{Input tree.} -\item{sim}{(list) : result of function \code{simulate}.} +\item{Y_data_vec}{: vector indicating the data at the tips.} -\item{Y_data}{: vector indicating the data at the tips.} +\item{sim}{(list) : result of function \code{simulate}.} -\item{Sigma_YY_inv}{: invert of the variance-covariance matrix of the data.} +\item{Sigma_YY_chol_inv}{: invert of the Cholesky variance-covariance matrix of the data.} } \description{ \code{compute_residuals.simple} computes the residuals after the fit of the diff --git a/man/enlight.Rd b/man/enlight.Rd index 635e5cc..6b56c32 100644 --- a/man/enlight.Rd +++ b/man/enlight.Rd @@ -30,9 +30,9 @@ call of \code{\link{PhyloEM}}. } \section{Methods (by class)}{ \itemize{ -\item \code{PhyloEM}: \code{\link{PhyloEM}} object -}} +\item \code{enlight(PhyloEM)}: \code{\link{PhyloEM}} object +}} \seealso{ \code{\link{PhyloEM}}, \code{\link{imputed_traits.PhyloEM}}, \code{\link{plot.PhyloEM}} diff --git a/man/equivalent_shifts_values.Rd b/man/equivalent_shifts_values.Rd index a4af15d..b58ed8b 100644 --- a/man/equivalent_shifts_values.Rd +++ b/man/equivalent_shifts_values.Rd @@ -15,10 +15,6 @@ equivalent_shifts_values(phylo, eq_shifts_edges, T_tree_ac, vect_Y, p) \item{T_tree_ac}{matrix of incidence of the tree, result of function \code{\link{incidence.matrix}}, actualized with coefficients computed by function \code{\link{incidence_matrix_actualization_factors}}.} - -\item{shifts}{a list of positions and values of original shifts.} - -\item{beta_0}{value of the original optimal value at the root.} } \value{ Named list, with "shifts_values" a matrix of shifts values diff --git a/man/estimateEM.Rd b/man/estimateEM.Rd index 332856d..5f9a240 100644 --- a/man/estimateEM.Rd +++ b/man/estimateEM.Rd @@ -10,9 +10,9 @@ estimateEM( Y_data_imp = Y_data, process = c("BM", "OU", "scOU", "rBM"), independent = FALSE, - tol_EM = list(variance = 10^(-2), value.root = 10^(-2), exp.root = 10^(-2), var.root - = 10^(-2), selection.strength = 10^(-2), normalized_half_life = 10^(-2), - log_likelihood = 10^(-2)), + tol_EM = list(variance = 10^(-2), value.root = 10^(-2), exp.root = 10^(-2), var.root = + 10^(-2), selection.strength = 10^(-2), normalized_half_life = 10^(-2), log_likelihood + = 10^(-2)), Nbr_It_Max = 500, method.variance = c("simple", "upward_downward"), method.init = c("default", "lasso"), @@ -28,10 +28,10 @@ estimateEM( max_selection.strength = 100, use_sigma_for_lasso = TRUE, max_triplet_number = 10000, - min_params = list(variance = 0, value.root = -10^(5), exp.root = -10^(5), var.root = - 0, selection.strength = 0), - max_params = list(variance = 10^(5), value.root = 10^(5), exp.root = 10^(5), var.root - = 10^(5), selection.strength = 10^(5)), + min_params = list(variance = 0, value.root = -10^(5), exp.root = -10^(5), var.root = 0, + selection.strength = 0), + max_params = list(variance = 10^(5), value.root = 10^(5), exp.root = 10^(5), var.root = + 10^(5), selection.strength = 10^(5)), var.init.root = diag(1, nrow(Y_data)), variance.init = diag(1, nrow(Y_data), nrow(Y_data)), methods.segmentation = c("lasso", "same_shifts", "best_single_move"), diff --git a/man/extract.variance_covariance.Rd b/man/extract.variance_covariance.Rd index 22e464b..88e84f2 100644 --- a/man/extract.variance_covariance.Rd +++ b/man/extract.variance_covariance.Rd @@ -15,12 +15,12 @@ \item{struct}{structural matrix of size (ntaxa+Nnode)*p, result of function \code{compute_variance_covariance}} -\item{what:}{sub-matrix to be extracted: +\item{what}{sub-matrix to be extracted: "YY" : sub-matrix of tips (p*ntaxa first lines and columns) "YZ" : sub matrix tips x nodes (p*Nnode last rows and p*ntaxa first columns) "ZZ" : sub matrix of nodes (p*Nnode last rows and columns)} -\item{miss;}{missing values of Y_data} +\item{miss}{missing values of Y_data} } \value{ sub-matrix of variance covariance. diff --git a/man/find_shift_values.Rd b/man/find_shift_values.Rd index 2954a62..43e4289 100644 --- a/man/find_shift_values.Rd +++ b/man/find_shift_values.Rd @@ -11,8 +11,6 @@ find_shift_values(shifts_edges, T_tree_ac, vect_Y, p, ntaxa) \item{T_tree_ac}{matrix of incidence of the tree, result of function \code{incidence.matrix}.} - -\item{m_Y}{the vector of values of the means at the tips.} } \value{ vector, with first entry the values at the root, and other entries the diff --git a/man/go_back_to_original_process.Rd b/man/go_back_to_original_process.Rd index 50563fe..3a65887 100644 --- a/man/go_back_to_original_process.Rd +++ b/man/go_back_to_original_process.Rd @@ -12,14 +12,14 @@ go_back_to_original_process( ) } \arguments{ -\item{phy_original:}{the original phylogenetic tree} +\item{phy_original}{the original phylogenetic tree} -\item{known.selection.strength:}{the known selection strength of the original +\item{known.selection.strength}{the known selection strength of the original OU.} -\item{sBM_variance:}{boolean. Is the root random ?} +\item{sBM_variance}{boolean. Is the root random ?} -\item{params:}{the inferred parameters of the BM on the re-scaled tree.} +\item{params}{the inferred parameters of the BM on the re-scaled tree.} } \value{ params_scOU the equivalent parameters of the OU on the original tree. diff --git a/man/imputed_traits.Rd b/man/imputed_traits.Rd index 6b28c84..c45cd40 100644 --- a/man/imputed_traits.Rd +++ b/man/imputed_traits.Rd @@ -58,9 +58,9 @@ state reconstruction) or at the tips (data imputation) } \section{Methods (by class)}{ \itemize{ -\item \code{PhyloEM}: \code{\link{PhyloEM}} object -}} +\item \code{imputed_traits(PhyloEM)}: \code{\link{PhyloEM}} object +}} \seealso{ \code{\link{params_process.PhyloEM}}, \code{\link{PhyloEM}} } diff --git a/man/init.parsimonyCost.Rd b/man/init.parsimonyCost.Rd index 1e7826a..332fe4b 100644 --- a/man/init.parsimonyCost.Rd +++ b/man/init.parsimonyCost.Rd @@ -7,9 +7,9 @@ init.parsimonyCost(phy, clusters) } \arguments{ -\item{clusters}{the vector of the clusters of the tips.} +\item{phy}{a phylogenetic tree, class \code{\link[ape]{phylo}}.} -\item{phylo}{a phylogenetic tree, class \code{\link[ape]{phylo}}.} +\item{clusters}{the vector of the clusters of the tips.} } \value{ A (ntaxa + Nnode)x(nclus) matrix, with ntaxa first lines initialized as diff --git a/man/init.simulate.OU.Rd b/man/init.simulate.OU.Rd index 4308439..050e51c 100644 --- a/man/init.simulate.OU.Rd +++ b/man/init.simulate.OU.Rd @@ -7,15 +7,15 @@ init.simulate.OU(phy, p, root.state, optimal.value, simulate_random, ...) } \arguments{ +\item{phy}{Input tree.} + +\item{p}{dimension of the trait simulated} + \item{root.state}{(list): state of the root, with: random : random state (TRUE) or deterministic state (FALSE) value.root : if deterministic, value of the character at the root exp.root : if random, expectation of the character at the root var.root : if random, variance of the character at the root (pxp matrix)} - -\item{phy:}{Input tree.} - -\item{p:}{dimension of the trait simulated} } \value{ paramSimu: array p x Nnode x 3, filled with NAs. diff --git a/man/init.simulate.StateAndExp.Rd b/man/init.simulate.StateAndExp.Rd index 719a319..5b4c18c 100644 --- a/man/init.simulate.StateAndExp.Rd +++ b/man/init.simulate.StateAndExp.Rd @@ -7,15 +7,15 @@ init.simulate.StateAndExp(phy, p, root.state, simulate_random) } \arguments{ +\item{phy}{Input tree.} + +\item{p}{dimension of the trait simulated} + \item{root.state}{(list): state of the root, with: random : random state (TRUE) or deterministic state (FALSE) value.root : if deterministic, value of the character at the root exp.root : if random, expectation of the character at the root var.root : if random, variance of the character at the root (pxp matrix)} - -\item{phy:}{Input tree.} - -\item{p:}{dimension of the trait simulated} } \value{ paramSimu: array p x Nnode x 2 (BM), filled with NAs. diff --git a/man/is.in.ranges.Rd b/man/is.in.ranges.Rd index f0464ae..d994a38 100644 --- a/man/is.in.ranges.Rd +++ b/man/is.in.ranges.Rd @@ -7,11 +7,11 @@ is.in.ranges(p, min, max) } \arguments{ +\item{p}{list of parameters with the correct structure} + \item{min}{list of minimum values for the parameters} \item{max}{list of maximum values for the parameters} - -\item{params}{list of parameters with the correct structure} } \value{ boolean diff --git a/man/log_likelihood.Rd b/man/log_likelihood.Rd index 8e76e2f..0d5ef46 100644 --- a/man/log_likelihood.Rd +++ b/man/log_likelihood.Rd @@ -34,11 +34,11 @@ The log likelihood of the data with the provided parameters on the tree. } \section{Methods (by class)}{ \itemize{ -\item \code{params_process}: \code{\link{params_process}} object +\item \code{log_likelihood(params_process)}: \code{\link{params_process}} object -\item \code{PhyloEM}: \code{\link{PhyloEM}} object -}} +\item \code{log_likelihood(PhyloEM)}: \code{\link{PhyloEM}} object +}} \seealso{ \code{\link{params_process}}, \code{\link{PhyloEM}} } diff --git a/man/merge_params_independent.Rd b/man/merge_params_independent.Rd index 14d6164..1d32ebf 100644 --- a/man/merge_params_independent.Rd +++ b/man/merge_params_independent.Rd @@ -7,7 +7,7 @@ merge_params_independent(params_split) } \arguments{ -\item{params_split:}{a list of parameters} +\item{params_split}{a list of parameters} } \value{ A parameter object diff --git a/man/model_selection.Rd b/man/model_selection.Rd index 8b486e0..a90e531 100644 --- a/man/model_selection.Rd +++ b/man/model_selection.Rd @@ -50,9 +50,9 @@ object, and returns the same fitted object. } \section{Methods (by class)}{ \itemize{ -\item \code{PhyloEM}: \code{\link{PhyloEM}} object -}} +\item \code{model_selection(PhyloEM)}: \code{\link{PhyloEM}} object +}} \seealso{ \code{\link{PhyloEM}}, \code{\link{params_process.PhyloEM}}, \code{\link{imputed_traits.PhyloEM}} diff --git a/man/penalty_pBIC.Rd b/man/penalty_pBIC.Rd index 82f98f6..2242d5e 100644 --- a/man/penalty_pBIC.Rd +++ b/man/penalty_pBIC.Rd @@ -24,9 +24,6 @@ penalty_pBIC( \item{K}{the dimension of the model.} \item{ntaxa}{the number of tips.} - -\item{C}{a constant, C > 1. Default is C = 1.1 -(as suggested in Baraud Giraud Huet (2009))} } \value{ value of the penalty. diff --git a/man/qr.solve_exact.Rd b/man/qr_solve_exact.Rd similarity index 89% rename from man/qr.solve_exact.Rd rename to man/qr_solve_exact.Rd index da90d0a..0484477 100644 --- a/man/qr.solve_exact.Rd +++ b/man/qr_solve_exact.Rd @@ -1,10 +1,10 @@ % Generated by roxygen2: do not edit by hand % Please edit documentation in R/parsimonyNumber.R -\name{qr.solve_exact} -\alias{qr.solve_exact} +\name{qr_solve_exact} +\alias{qr_solve_exact} \title{exact qr.solve} \usage{ -qr.solve_exact(a, b, tol = 1e-07) +qr_solve_exact(a, b, tol = 1e-07) } \arguments{ \item{a}{a QR decomposition or a rectangular matrix.} diff --git a/man/segmentation.OU.specialCase.lasso.Rd b/man/segmentation.OU.specialCase.lasso.Rd index f61470f..84a31bc 100644 --- a/man/segmentation.OU.specialCase.lasso.Rd +++ b/man/segmentation.OU.specialCase.lasso.Rd @@ -17,11 +17,6 @@ segmentation.OU.specialCase.lasso( \item{phylo}{a phylogenetic tree} \item{nbr_of_shifts}{Number of shifts on the phylogeny allowed} - -\item{conditional_law_X}{moments of the conditional law of X given Y, result -of function \code{compute_M.OU.specialCase}} - -\item{selection.strength}{the selection strength} } \value{ List containing : beta_0 : the optimal value at the root diff --git a/man/segmentation.OU.specialCase.same_shifts.Rd b/man/segmentation.OU.specialCase.same_shifts.Rd index e333791..ff72680 100644 --- a/man/segmentation.OU.specialCase.same_shifts.Rd +++ b/man/segmentation.OU.specialCase.same_shifts.Rd @@ -11,11 +11,6 @@ segmentation.OU.specialCase.same_shifts(phylo, shifts_old, D, Xp, ...) \item{phylo}{a phylogenetic tree} \item{shifts_old}{the previous list of shifts (only position is used)} - -\item{conditional_law_X}{moments of the conditional law of X given Y, result -of function \code{compute_M.OU.specialCase}} - -\item{selection.strength}{the selection strength} } \value{ List containing : beta_0 : the optimal value at the root diff --git a/man/simul_process.Rd b/man/simul_process.Rd index 353434d..da3b3a9 100644 --- a/man/simul_process.Rd +++ b/man/simul_process.Rd @@ -66,11 +66,11 @@ An S3 object of class \code{simul_process}. This contains: } \section{Methods (by class)}{ \itemize{ -\item \code{params_process}: \code{\link{params_process}} object +\item \code{simul_process(params_process)}: \code{\link{params_process}} object -\item \code{PhyloEM}: \code{\link{PhyloEM}} object -}} +\item \code{simul_process(PhyloEM)}: \code{\link{PhyloEM}} object +}} \seealso{ \code{\link{params_process}}, \code{\link{PhyloEM}}, \code{\link{extract.simul_process}} } diff --git a/man/simulate_internal.Rd b/man/simulate_internal.Rd index 9885818..9ef33fe 100644 --- a/man/simulate_internal.Rd +++ b/man/simulate_internal.Rd @@ -8,8 +8,8 @@ simulate_internal( phylo, process = c("BM", "OU", "scOU", "OUBM", "StudentOU"), p = 1, - root.state = list(random = FALSE, stationary.root = FALSE, value.root = NA, exp.root - = NA, var.root = NA), + root.state = list(random = FALSE, stationary.root = FALSE, value.root = NA, exp.root = + NA, var.root = NA), shifts = list(edges = NULL, values = NULL, relativeTimes = NULL), eps = 10^(-6), selection.strength = NULL, diff --git a/tests/testthat/test-utilities-rotations.R b/tests/testthat/test-utilities-rotations.R index d953868..a279a6d 100644 --- a/tests/testthat/test-utilities-rotations.R +++ b/tests/testthat/test-utilities-rotations.R @@ -131,10 +131,11 @@ test_that("Rotations", { # rotate params expect_equivalent(res$alpha_max$params_estim$`3`, rotate_params(res_rot$alpha_max$params_estim$`3`, rot), - 1e-5) + 1e-3) expect_equal(log_likelihood(params_process(res, K = 3), phylo = monkeys$phy, Y_data = monkeys$dat), - log_likelihood(rotate_params(params_process(res_rot, K = 3), rot), phylo = monkeys$phy, Y_data = monkeys$dat)) + log_likelihood(rotate_params(params_process(res_rot, K = 3), rot), phylo = monkeys$phy, Y_data = monkeys$dat), + tolerance = 1e-4) # another rotation From 2b7bf995887ac865d2a321bcd883ef9a552a474d Mon Sep 17 00:00:00 2001 From: Paul Bastide Date: Thu, 5 Oct 2023 10:17:44 +0200 Subject: [PATCH 03/10] node_optimal_values function --- NAMESPACE | 1 + R/shifts_manipulations.R | 63 ++++++++++++++++++++++++++++++-- man/compute_betas_from_shifts.Rd | 5 +++ man/node_optimal_values.Rd | 41 +++++++++++++++++++++ 4 files changed, 107 insertions(+), 3 deletions(-) create mode 100644 man/node_optimal_values.Rd diff --git a/NAMESPACE b/NAMESPACE index aa3f0d8..8471ef9 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -53,6 +53,7 @@ export(incidence.matrix.full) export(log_likelihood) export(merge_rotations) export(model_selection) +export(node_optimal_values) export(params_BM) export(params_OU) export(params_process) diff --git a/R/shifts_manipulations.R b/R/shifts_manipulations.R index f866d82..ddbbfc1 100644 --- a/R/shifts_manipulations.R +++ b/R/shifts_manipulations.R @@ -474,6 +474,12 @@ update.compute_betas_from_shifts <- function(edgeNbr, ancestral, shifts, ...){ #' \code{compute_betas_from_shifts} computes the optimal values at the nodes and tips of the #' tree, given the value at the root and the list of shifts occurring in the tree. #' It assumes an OU model. +#' +#' @details +#' Note that this is intended to be an internal function, and should not be used. +#' In general, use \code{\link{node_optimal_values}} to get optimal values +#' from a set of parameters. +#' #' # @details # This function uses function \code{recursionDown} for recursion on the tree, @@ -487,9 +493,6 @@ update.compute_betas_from_shifts <- function(edgeNbr, ancestral, shifts, ...){ #' @return Vector of size (ntaxa + Nnode) of the optimal values at the tips #' of the tree. #' -#' @example -#' -#' #' @export #' #06/10/14 - Initial release @@ -504,6 +507,60 @@ compute_betas_from_shifts <- function(phylo, optimal.value, shifts){ return(betas) } +## +#' @title Computation of the optimal values at nodes and tips. +#' +#' @description +#' \code{compute_betas_from_shifts} computes the optimal values at the nodes and tips of the +#' tree, given the value at the root and the list of shifts occurring in the tree. +#' It assumes an OU model. +# +#' @param phylo a phylogenetic tree, class \code{\link[ape]{phylo}}. +#' @param param an object of class \code{\link{params_process}}. +#' +#' @return Matrix of size ntraits x (ntaxa + Nnode) of the optimal values at +#' the node and tips of the tree. +#' Column names correspond to the number of the node in the phylo object. +#' +#' @examples +#' set.seed(1792) +#' ntaxa = 10 +#' tree <- rphylo(ntaxa, 1, 0.1) +#' # parameters of the process +#' par <- params_process("BM", ## Process +#' p = 2, ## Dimension +#' variance = diag(0.5, 2, 2) + 0.5, ## Rate matrix +#' edges = c(4, 10, 15), ## Positions of the shifts +#' values = cbind(c(5, 4), ## Values of the shifts +#' c(-4, -5), +#' c(5, -3))) +#' plot(par, phylo = tree, traits = 1, value_in_box = TRUE, +#' shifts_bg = "white", root_bg = "white", ancestral_as_shift = TRUE, root_adj = 5) +#' nodelabels() +#' node_optimal_values(par, tree) +#' +#' @export +#' +## +node_optimal_values <- function(param, phylo){ + if(param$root.state$random) { + opt_root <- param$root.state$exp.root + } else { + opt_root <- param$root.state$value.root + } + p <- ncol(param$variance) + if (is.null(p)) p <- 1 + res <- matrix(NA, nrow = p, ncol = Nnode(phylo) + Ntip(phylo)) + for (dim in 1:ncol(param$variance)) { + shifts <- param$shifts + shifts$values <- shifts$values[dim, ] + res[dim, ] <- compute_betas_from_shifts(phylo, opt_root[dim], shifts) + } + colnames(res) <- 1:ncol(res) + rownames(res) <- colnames(param$variance) + return(res) +} + ## #' @title Initialization for the allocation of shifts. #' diff --git a/man/compute_betas_from_shifts.Rd b/man/compute_betas_from_shifts.Rd index 4f76ccc..fbc6609 100644 --- a/man/compute_betas_from_shifts.Rd +++ b/man/compute_betas_from_shifts.Rd @@ -22,3 +22,8 @@ of the tree. tree, given the value at the root and the list of shifts occurring in the tree. It assumes an OU model. } +\details{ +Note that this is intended to be an internal function, and should not be used. +In general, use \code{\link{node_optimal_values}} to get optimal values +from a set of parameters. +} diff --git a/man/node_optimal_values.Rd b/man/node_optimal_values.Rd new file mode 100644 index 0000000..4a4fc59 --- /dev/null +++ b/man/node_optimal_values.Rd @@ -0,0 +1,41 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/shifts_manipulations.R +\name{node_optimal_values} +\alias{node_optimal_values} +\title{Computation of the optimal values at nodes and tips.} +\usage{ +node_optimal_values(param, phylo) +} +\arguments{ +\item{param}{an object of class \code{\link{params_process}}.} + +\item{phylo}{a phylogenetic tree, class \code{\link[ape]{phylo}}.} +} +\value{ +Matrix of size ntraits x (ntaxa + Nnode) of the optimal values at +the node and tips of the tree. +Column names correspond to the number of the node in the phylo object. +} +\description{ +\code{compute_betas_from_shifts} computes the optimal values at the nodes and tips of the +tree, given the value at the root and the list of shifts occurring in the tree. +It assumes an OU model. +} +\examples{ +set.seed(1792) +ntaxa = 10 +tree <- rphylo(ntaxa, 1, 0.1) +# parameters of the process +par <- params_process("BM", ## Process + p = 2, ## Dimension + variance = diag(0.5, 2, 2) + 0.5, ## Rate matrix + edges = c(4, 10, 15), ## Positions of the shifts + values = cbind(c(5, 4), ## Values of the shifts + c(-4, -5), + c(5, -3))) +plot(par, phylo = tree, traits = 1, value_in_box = TRUE, + shifts_bg = "white", root_bg = "white", ancestral_as_shift = TRUE, root_adj = 5) +nodelabels() +node_optimal_values(par, tree) + +} From 61d007a1eebcc8f68893761f1d92312f17885a91 Mon Sep 17 00:00:00 2001 From: Paul Bastide Date: Fri, 6 Oct 2023 13:00:04 +0200 Subject: [PATCH 04/10] Adding trait names in params_process. --- R/estimateEM.R | 32 ++++++++++++++++++++++++++++++++ R/init_EM.R | 12 ++++++++++++ man/params_BM.Rd | 3 +++ man/params_OU.Rd | 3 +++ 4 files changed, 50 insertions(+) diff --git a/R/estimateEM.R b/R/estimateEM.R index 163588f..b930f97 100644 --- a/R/estimateEM.R +++ b/R/estimateEM.R @@ -1425,6 +1425,7 @@ params_process.PhyloEM <- function(x, method.selection = NULL, variance = res$variance, selection.strength = res$selection.strength, optimal.value = res$optimal.value) + if (!is.null(rownames(x$Y_data))) res <- name_params(res, rownames(x$Y_data)) res$variance <- as(res$variance, "dpoMatrix") class(res) <- "params_process" if (attr(res, "Neq") > 1){ @@ -1433,6 +1434,37 @@ params_process.PhyloEM <- function(x, method.selection = NULL, return(res) } +name_params <- function(res, names) { + ## root state + if (isnonnullna(res$root.state$value.root)) names(res$root.state$value.root) <- names + if (isnonnullna(res$root.state$exp.root)) names(res$root.state$exp.root) <- names + res$root.state$var.root <- name_matrix(res$root.state$var.root, names) + ## shifts + if (isnonnullna(res$shifts$values)) rownames(res$shifts$values) <- names + ## variance + res$variance <- name_matrix(res$variance, names) + ## selection strength + res$selection.strength <- name_matrix(res$selection.strength, names) + ## optimal values + if (isnonnullna(res$optimal.value)) names(res$optimal.value) <- names + return(res) +} + +isnonnullna <- function(x) { + return(!is.null(x) && !any(is.na(x))) +} + +name_matrix <- function(M, names) { + if (isnonnullna(M)) { + if (!is.null(attr(class(M), "package")) && attr(class(M), "package") == "Matrix") { + dimnames(M) <- rep.int(list(names), 2L) + } else { + colnames(M) <- rownames(M) <- names + } + } + return(M) +} + extract_params <- function(x, method, alpha_str){ if (!is.null(x[[alpha_str]][[method]])){ res <- x[[alpha_str]][[method]]$params_select diff --git a/R/init_EM.R b/R/init_EM.R index f68ed67..38549aa 100644 --- a/R/init_EM.R +++ b/R/init_EM.R @@ -351,6 +351,7 @@ params_process.character <- function(x, ...){ #' @param sBM_variance if the root is random, does it depend on the length of the #' root edge ? (For equivalent purposes with a rescaled OU). Default to FALSE. If #' TRUE, a phylogenetic tree with root edge length must be provided. +#' @param trait_names vector of trait names values. Must be of length p. #' @param ... unused. #' #' @return an object of class \code{params_process}. @@ -372,6 +373,7 @@ params_BM <- function(p = 1, nbr_of_shifts = length(edges), phylo = NULL, sBM_variance = FALSE, + trait_names = NULL, ...) { if (random) { value.root <- NA @@ -415,6 +417,10 @@ params_BM <- function(p = 1, params_init$root.state <- test.root.state(params_init$root.state, "BM") params_init$variance <- as(params_init$variance, "dpoMatrix") params_init$process <- "BM" + if (!is.null(trait_names)) { + if(length(trait_names) != p) stop("`trait_names` must be a vector of dimension p.") + params_init <- name_params(params_init, trait_names) + } class(params_init) <- "params_process" return(params_init) } @@ -459,6 +465,7 @@ params_BM <- function(p = 1, #' provided (to allow a random sampling of its edges). #' @param phylo a phylogenetic tree of class \code{\link[ape]{phylo}}. Needed only if #' the shifts edges are not specified. +#' @param trait_names vector of trait names values. Must be of length p. #' @param ... unused. #' #' @return an object of class \code{params_process}. @@ -482,6 +489,7 @@ params_OU <- function(p = 1, relativeTimes = NULL, nbr_of_shifts = length(edges), phylo = NULL, + trait_names = NULL, ...) { if (random) { value.root <- NA @@ -533,6 +541,10 @@ params_OU <- function(p = 1, selection.strength = params_init$selection.strength, optimal.value = optimal.value) params_init$variance <- as(params_init$variance, "dpoMatrix") + if (!is.null(trait_names)) { + if(length(trait_names) != p) stop("`trait_names` must be a vector of dimension p.") + params_init <- name_params(params_init, trait_names) + } params_init$process <- "OU" class(params_init) <- "params_process" return(params_init) diff --git a/man/params_BM.Rd b/man/params_BM.Rd index 110174c..d71cbc3 100644 --- a/man/params_BM.Rd +++ b/man/params_BM.Rd @@ -17,6 +17,7 @@ params_BM( nbr_of_shifts = length(edges), phylo = NULL, sBM_variance = FALSE, + trait_names = NULL, ... ) } @@ -60,6 +61,8 @@ length (slot root.edge).} root edge ? (For equivalent purposes with a rescaled OU). Default to FALSE. If TRUE, a phylogenetic tree with root edge length must be provided.} +\item{trait_names}{vector of trait names values. Must be of length p.} + \item{...}{unused.} } \value{ diff --git a/man/params_OU.Rd b/man/params_OU.Rd index 14bf1f6..824d88d 100644 --- a/man/params_OU.Rd +++ b/man/params_OU.Rd @@ -19,6 +19,7 @@ params_OU( relativeTimes = NULL, nbr_of_shifts = length(edges), phylo = NULL, + trait_names = NULL, ... ) } @@ -69,6 +70,8 @@ provided (to allow a random sampling of its edges).} \item{phylo}{a phylogenetic tree of class \code{\link[ape]{phylo}}. Needed only if the shifts edges are not specified.} +\item{trait_names}{vector of trait names values. Must be of length p.} + \item{...}{unused.} } \value{ From 59856835ddc897471d1c1d059b27e0ebdb33923e Mon Sep 17 00:00:00 2001 From: Paul Bastide Date: Fri, 6 Oct 2023 13:55:09 +0200 Subject: [PATCH 05/10] Update GHA --- .github/workflows/R-CMD-check-no-force-suggest.yaml | 6 +++--- .github/workflows/check-standard.yaml | 10 +++++----- .github/workflows/pkgdown.yaml | 8 ++++---- .github/workflows/test-coverage.yaml | 4 ++-- 4 files changed, 14 insertions(+), 14 deletions(-) diff --git a/.github/workflows/R-CMD-check-no-force-suggest.yaml b/.github/workflows/R-CMD-check-no-force-suggest.yaml index 83e8e64..03c6228 100644 --- a/.github/workflows/R-CMD-check-no-force-suggest.yaml +++ b/.github/workflows/R-CMD-check-no-force-suggest.yaml @@ -27,11 +27,11 @@ jobs: steps: - uses: actions/checkout@v2 - - uses: r-lib/actions/setup-r@master + - uses: r-lib/actions/setup-r@v2 with: r-version: ${{ matrix.config.r }} - - uses: r-lib/actions/setup-pandoc@master + - uses: r-lib/actions/setup-pandoc@v2 - name: Query dependencies run: | @@ -72,7 +72,7 @@ jobs: - name: Upload check results if: failure() - uses: actions/upload-artifact@master + uses: actions/upload-artifact@v3 with: name: ${{ runner.os }}-r${{ matrix.config.r }}-results path: check diff --git a/.github/workflows/check-standard.yaml b/.github/workflows/check-standard.yaml index 4d2915d..cd38fd1 100644 --- a/.github/workflows/check-standard.yaml +++ b/.github/workflows/check-standard.yaml @@ -29,21 +29,21 @@ jobs: R_KEEP_PKG_SOURCE: yes steps: - - uses: actions/checkout@v2 + - uses: actions/checkout@v3 - - uses: r-lib/actions/setup-pandoc@v1 + - uses: r-lib/actions/setup-pandoc@v2 - - uses: r-lib/actions/setup-r@v1 + - uses: r-lib/actions/setup-r@v2 with: r-version: ${{ matrix.config.r }} http-user-agent: ${{ matrix.config.http-user-agent }} use-public-rspm: true - - uses: r-lib/actions/setup-r-dependencies@v1 + - uses: r-lib/actions/setup-r-dependencies@v2 with: extra-packages: rcmdcheck - - uses: r-lib/actions/check-r-package@v1 + - uses: r-lib/actions/check-r-package@v2 - name: Show testthat output if: always() diff --git a/.github/workflows/pkgdown.yaml b/.github/workflows/pkgdown.yaml index 59ae308..65b08c8 100644 --- a/.github/workflows/pkgdown.yaml +++ b/.github/workflows/pkgdown.yaml @@ -13,15 +13,15 @@ jobs: env: GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} steps: - - uses: actions/checkout@v2 + - uses: actions/checkout@v3 - - uses: r-lib/actions/setup-pandoc@v1 + - uses: r-lib/actions/setup-pandoc@v2 - - uses: r-lib/actions/setup-r@v1 + - uses: r-lib/actions/setup-r@v2 with: use-public-rspm: true - - uses: r-lib/actions/setup-r-dependencies@v1 + - uses: r-lib/actions/setup-r-dependencies@v2 with: extra-packages: pkgdown needs: website diff --git a/.github/workflows/test-coverage.yaml b/.github/workflows/test-coverage.yaml index 084a108..cab4154 100644 --- a/.github/workflows/test-coverage.yaml +++ b/.github/workflows/test-coverage.yaml @@ -15,9 +15,9 @@ jobs: GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} steps: - - uses: actions/checkout@v2 + - uses: actions/checkout@v3 - - uses: r-lib/actions/setup-r@v1 + - uses: r-lib/actions/setup-r@v2 with: use-public-rspm: true From 0493bf81207f803bcd1e81f98b4c46c780a17186 Mon Sep 17 00:00:00 2001 From: Paul Bastide Date: Wed, 31 Jan 2024 14:22:54 +0100 Subject: [PATCH 06/10] Documentation issues --- NEWS.md | 5 ++ R/M_step.R | 2 +- R/generic_functions.R | 3 +- R/parsimonyNumber.R | 4 +- R/shifts_manipulations.R | 6 +- R/shutoff.R | 80 ------------------- ...ompute_actualization_matrix_ultrametric.Rd | 2 +- man/compute_sum_var_diff.Rd | 2 +- man/find_grid_alpha.Rd | 3 +- man/incidence_matrix_actualization_factors.Rd | 2 +- man/init.compute_betas_from_shifts.Rd | 2 +- man/init.parsimonyCost.Rd | 2 +- man/init.parsimonyNumber.Rd | 2 +- man/is.in.ranges.Rd | 63 --------------- 14 files changed, 21 insertions(+), 157 deletions(-) diff --git a/NEWS.md b/NEWS.md index b7f5cce..714e5d8 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,4 +1,9 @@ # PhylogeneticEM 1.6.0.9000 +* Bug fixes: + * Fix documentation notes on CRAN checks. +* New functions and features: + * Function `node_optimal_values` for optimal values computation at each nodes. + * Trait names are now inherited by the `params_process` objects. # PhylogeneticEM 1.6.0 * Deprecation fix: diff --git a/R/M_step.R b/R/M_step.R index 5d95d78..97c438a 100644 --- a/R/M_step.R +++ b/R/M_step.R @@ -543,7 +543,7 @@ compute_var_diff.OU <- function(phylo, conditional_law_X, selection.strength) { #' @title Compute weighted sum of var_diff #' #' @description -#' \code{compute_sum_var_diff} computes sum_{e edge} ell_j * Var[X_j - X_pa(j) | Y] +#' \code{compute_sum_var_diff} computes \sum_\{e edge\} ell_j * Var[X_j - X_pa(j) | Y] #' #' @param phylo a phylogenetic tree #' @param var_diff result of function \code{compute_var_diff.BM} diff --git a/R/generic_functions.R b/R/generic_functions.R index 8104782..85f0b35 100644 --- a/R/generic_functions.R +++ b/R/generic_functions.R @@ -666,7 +666,8 @@ check_dimensions.shifts <- function(p, shifts){ #' @title Find a reasonable grid for alpha #' #' @description Grid so that -#' 2*ln(2)*quantile(d_ij)/factor_up_alpha < t_{1/2} < factor_down_alpha * ln(2) * h_tree +#' 2*ln(2)*quantile(d_ij)/factor_up_alpha < t_1/2 < factor_down_alpha * ln(2) * h_tree, +#' with t_1/2 the phylogenetic half life: t_1/2 = log(2)/alpha. #' Ensures that for alpha_min, it is almost a BM, and for alpha_max, #' almost all the tips are decorrelated. #' diff --git a/R/parsimonyNumber.R b/R/parsimonyNumber.R index 9654981..9bf461e 100644 --- a/R/parsimonyNumber.R +++ b/R/parsimonyNumber.R @@ -150,7 +150,7 @@ extract.parsimonyCost <- function(x, #' #' @details #' At a tip i in state k, the line-vector is initialized as follow : -#' (1 - Ind(k=p)_{1<=p<=nclus})*Inf (where Inf * 0 = 0) +#' (1 - Ind(k=p)_\{1<=p<=nclus\})*Inf (where Inf * 0 = 0) #' #' @param phy a phylogenetic tree, class \code{\link[ape]{phylo}}. #' @param clusters the vector of the clusters of the tips. @@ -356,7 +356,7 @@ extract.parsimonyNumber <- function(x, #' NAs everywhere, except for the tips. #' #' @details -#' At a tip i in state k, the line-vector is initialized as follow : Ind(k=p)_{1<=p<=nclus} +#' At a tip i in state k, the line-vector is initialized as follow : Ind(k=p)_\{1<=p<=nclus\} #' #' @param phy phylogenetic tree. #' @param clusters the vector of the clusters of the tips. diff --git a/R/shifts_manipulations.R b/R/shifts_manipulations.R index ddbbfc1..1309a15 100644 --- a/R/shifts_manipulations.R +++ b/R/shifts_manipulations.R @@ -334,7 +334,7 @@ shifts.matrix_to_list <- function(delta){ #' #' @description #' \code{incidence_matrix_actualization_factors} computes a ntaxa x Nedge matrix of the -#' (1 - exp(-alpha * (t_i - t_pa(j) - nu_j * l_j)))_{i tip, j node}. +#' (1 - exp(-alpha * (t_i - t_pa(j) - nu_j * l_j)))_\{i tip, j node\}. #' This matrix is to be multiplied to the incidence matrix with an outer product. #' #' @param tree a phylogenetic tree. @@ -375,7 +375,7 @@ incidence_matrix_actualization_factors <- function(tree, #' #' @description #' \code{compute_actualization_matrix_ultrametric} computes a squares p*Nedge bloc diagonal -#' matrix of the (I_p - exp(-A * (h - t_pa(j))))_{j node}. +#' matrix of the (I_p - exp(-A * (h - t_pa(j))))_\{j node\}. #' #' @details #' Careful: the root is not taken into account in this function. @@ -419,7 +419,7 @@ compute_actualization_matrix_ultrametric <- function(tree, #' tips, with the value at the root. #' #' @details -#' This function is used in function \code{compute_betas_from_shifts_from_shifts} and is designed to +#' This function is used in function \code{compute_betas_from_shifts} and is designed to #' furnish function \code{update.compute_betas_from_shifts} with the right structure of data. #' #' @param phy Input tree. diff --git a/R/shutoff.R b/R/shutoff.R index 0050fe9..d6876fe 100644 --- a/R/shutoff.R +++ b/R/shutoff.R @@ -120,86 +120,6 @@ shutoff.EM.OU.stationary.root_AND_shifts_at_nodes.half_life <- function(params_o } } -#' ###################################################### -#' ## Finite parameters ? -#' ###################################################### -# ## -# # @title Check whether parameters are finite. -# # -# # @description -# # \code{is.finite.params} checks whether calculated parameters in the EM are -# # finite or not. -# # -# # @details -# # This function is used to test the convergence of the algorithm. -# # -# # @param params list of parameters with the correct structure -# # -# # @return boolean -# # -# # @keywords internal -# # -# #10/06/14 - Initial release -# ## -#' -#' is.finite.params.BM <- function(params) { -#' if (params$root.state$random) { -#' return(is.finite.params.BM.randroot(params)) -#' } else { -#' return(is.finite.params.BM.fixedroot(params)) -#' } -#' } -#' -#' is.finite.params.BM.randroot <- function(params) { -#' if (is.finite(params$variance) && -#' is.finite(params$root.state$exp.root) && -#' is.finite(params$root.state$var.root) ) { -#' return(TRUE) -#' } else { -#' return(FALSE) -#' } -#' } -#' -#' is.finite.params.BM.fixedroot <- function(params) { -#' if ( is.finite(params$variance) && -#' is.finite(params$root.state$value.root) ) { -#' return(TRUE) -#' } else { -#' return(FALSE) -#' } -#' } -#' -#' is.finite.params.OU <- function(stationary.root, shifts_at_nodes, alpha_known){ -#' if (stationary.root && shifts_at_nodes && alpha_known) { -#' return(is.finite.params.OU.specialCase) -#' } else if (stationary.root && shifts_at_nodes) { -#' return(is.finite.params.OU.stationary.root_AND_shifts_at_nodes) -#' } else { -#' stop("The EM algorithm for the OU is only defined (for the moment) for a stationary root and shifts at nodes !") -#' } -#' } -#' -#' is.finite.params.OU.specialCase <- function(params) { -#' if (is.finite(params$variance) && -#' is.finite(params$root.state$exp.root) && -#' is.finite(params$root.state$var.root) ) { -#' return(TRUE) -#' } else { -#' return(FALSE) -#' } -#' } -#' -#' is.finite.params.OU.stationary.root_AND_shifts_at_nodes <- function(params) { -#' if (is.finite(params$variance) && -#' is.finite(params$root.state$exp.root) && -#' is.finite(params$root.state$var.root) && -#' is.finite(params$selection.strength)) { -#' return(TRUE) -#' } else { -#' return(FALSE) -#' } -#' } - ################################################################# ## Are parameters in ranges ? ################################################################# diff --git a/man/compute_actualization_matrix_ultrametric.Rd b/man/compute_actualization_matrix_ultrametric.Rd index 4ec1e57..09da6ca 100644 --- a/man/compute_actualization_matrix_ultrametric.Rd +++ b/man/compute_actualization_matrix_ultrametric.Rd @@ -22,7 +22,7 @@ Matrix of size p*Nedge } \description{ \code{compute_actualization_matrix_ultrametric} computes a squares p*Nedge bloc diagonal -matrix of the (I_p - exp(-A * (h - t_pa(j))))_{j node}. +matrix of the (I_p - exp(-A * (h - t_pa(j))))_\{j node\}. } \details{ Careful: the root is not taken into account in this function. diff --git a/man/compute_sum_var_diff.Rd b/man/compute_sum_var_diff.Rd index 4d51202..e2bd43a 100644 --- a/man/compute_sum_var_diff.Rd +++ b/man/compute_sum_var_diff.Rd @@ -15,6 +15,6 @@ compute_sum_var_diff(phylo, var_diff) matrix p x p } \description{ -\code{compute_sum_var_diff} computes sum_{e edge} ell_j * Var[X_j - X_pa(j) | Y] +\code{compute_sum_var_diff} computes \sum_\{e edge\} ell_j * Var[X_j - X_pa(j) | Y] } \keyword{internal} diff --git a/man/find_grid_alpha.Rd b/man/find_grid_alpha.Rd index a12064f..d33d6ba 100644 --- a/man/find_grid_alpha.Rd +++ b/man/find_grid_alpha.Rd @@ -42,7 +42,8 @@ A grid of alpha values } \description{ Grid so that -2*ln(2)*quantile(d_ij)/factor_up_alpha < t_{1/2} < factor_down_alpha * ln(2) * h_tree +2*ln(2)*quantile(d_ij)/factor_up_alpha < t_1/2 < factor_down_alpha * ln(2) * h_tree, +with t_1/2 the phylogenetic half life: t_1/2 = log(2)/alpha. Ensures that for alpha_min, it is almost a BM, and for alpha_max, almost all the tips are decorrelated. } diff --git a/man/incidence_matrix_actualization_factors.Rd b/man/incidence_matrix_actualization_factors.Rd index 1f9c1bf..599ffea 100644 --- a/man/incidence_matrix_actualization_factors.Rd +++ b/man/incidence_matrix_actualization_factors.Rd @@ -25,7 +25,7 @@ Matrix of size ntaxa x Nedge } \description{ \code{incidence_matrix_actualization_factors} computes a ntaxa x Nedge matrix of the -(1 - exp(-alpha * (t_i - t_pa(j) - nu_j * l_j)))_{i tip, j node}. +(1 - exp(-alpha * (t_i - t_pa(j) - nu_j * l_j)))_\{i tip, j node\}. This matrix is to be multiplied to the incidence matrix with an outer product. } \keyword{internal} diff --git a/man/init.compute_betas_from_shifts.Rd b/man/init.compute_betas_from_shifts.Rd index de22c24..e106320 100644 --- a/man/init.compute_betas_from_shifts.Rd +++ b/man/init.compute_betas_from_shifts.Rd @@ -20,7 +20,7 @@ Matrix of size (Nnode + ntaxa)x1 of NAs, with the optimal value tips, with the value at the root. } \details{ -This function is used in function \code{compute_betas_from_shifts_from_shifts} and is designed to +This function is used in function \code{compute_betas_from_shifts} and is designed to furnish function \code{update.compute_betas_from_shifts} with the right structure of data. } \keyword{internal} diff --git a/man/init.parsimonyCost.Rd b/man/init.parsimonyCost.Rd index 332fe4b..1de3959 100644 --- a/man/init.parsimonyCost.Rd +++ b/man/init.parsimonyCost.Rd @@ -21,6 +21,6 @@ NAs everywhere, except for the tips. } \details{ At a tip i in state k, the line-vector is initialized as follow : -(1 - Ind(k=p)_{1<=p<=nclus})*Inf (where Inf * 0 = 0) +(1 - Ind(k=p)_\{1<=p<=nclus\})*Inf (where Inf * 0 = 0) } \keyword{internal} diff --git a/man/init.parsimonyNumber.Rd b/man/init.parsimonyNumber.Rd index b6ae1fc..9c5f908 100644 --- a/man/init.parsimonyNumber.Rd +++ b/man/init.parsimonyNumber.Rd @@ -20,6 +20,6 @@ as described. NAs everywhere, except for the tips. } \details{ -At a tip i in state k, the line-vector is initialized as follow : Ind(k=p)_{1<=p<=nclus} +At a tip i in state k, the line-vector is initialized as follow : Ind(k=p)_\{1<=p<=nclus\} } \keyword{internal} diff --git a/man/is.in.ranges.Rd b/man/is.in.ranges.Rd index d994a38..bc6b54b 100644 --- a/man/is.in.ranges.Rd +++ b/man/is.in.ranges.Rd @@ -21,69 +21,6 @@ boolean in the defined ranges. } \details{ -###################################################### -## Finite parameters ? -###################################################### - -is.finite.params.BM <- function(params) { - if (params$root.state$random) { - return(is.finite.params.BM.randroot(params)) - } else { - return(is.finite.params.BM.fixedroot(params)) - } -} - -is.finite.params.BM.randroot <- function(params) { - if (is.finite(params$variance) && - is.finite(params$root.state$exp.root) && - is.finite(params$root.state$var.root) ) { - return(TRUE) - } else { - return(FALSE) - } -} - -is.finite.params.BM.fixedroot <- function(params) { - if ( is.finite(params$variance) && - is.finite(params$root.state$value.root) ) { - return(TRUE) - } else { - return(FALSE) - } -} - -is.finite.params.OU <- function(stationary.root, shifts_at_nodes, alpha_known){ - if (stationary.root && shifts_at_nodes && alpha_known) { - return(is.finite.params.OU.specialCase) - } else if (stationary.root && shifts_at_nodes) { - return(is.finite.params.OU.stationary.root_AND_shifts_at_nodes) - } else { - stop("The EM algorithm for the OU is only defined (for the moment) for a stationary root and shifts at nodes !") - } -} - -is.finite.params.OU.specialCase <- function(params) { - if (is.finite(params$variance) && - is.finite(params$root.state$exp.root) && - is.finite(params$root.state$var.root) ) { - return(TRUE) - } else { - return(FALSE) - } -} - -is.finite.params.OU.stationary.root_AND_shifts_at_nodes <- function(params) { - if (is.finite(params$variance) && - is.finite(params$root.state$exp.root) && - is.finite(params$root.state$var.root) && - is.finite(params$selection.strength)) { - return(TRUE) - } else { - return(FALSE) - } -} - - This function is used to test the convergence of the algorithm. } \keyword{internal} From b4e150eaeba6902bb984b752c8817004001ec86d Mon Sep 17 00:00:00 2001 From: Paul Bastide Date: Wed, 31 Jan 2024 15:46:24 +0100 Subject: [PATCH 07/10] Rename internal functions to avoid ambiguities. --- R/E_step.R | 40 ++++++++--------- R/M_step.R | 2 +- R/init_EM.R | 2 +- R/model_selection.R | 6 +-- R/partitionsNumber.R | 44 +++++++++---------- R/plot_functions.R | 4 +- R/simulate.R | 2 +- man/compute_E.simple.Rd | 2 +- man/compute_log_likelihood.simple.Rd | 2 +- man/compute_sum_var_diff.Rd | 2 +- ...ance.Rd => extract_variance_covariance.Rd} | 10 ++--- man/{prod.index.Rd => prod_index.Rd} | 10 ++--- man/{sum.partitions.Rd => sum_partitions.Rd} | 10 ++--- ...rdFixed.Rd => sum_partitions.cardFixed.Rd} | 10 ++--- man/{sum.prod.comb.Rd => sum_prod.comb.Rd} | 10 ++--- man/{sum.simplex.Rd => sum_simplex.Rd} | 10 ++--- man/update.partitionsNumber.gen.Rd | 2 +- .../multivariate_data_generation.R | 26 +++++------ .../multivariate_data_generation_two_bases.R | 26 +++++------ ...te_data_generation_two_bases_with errors.R | 26 +++++------ simulations_study/multivariate_exploitation.R | 8 ++-- simulations_study/penalty_calibration.R | 2 +- simulations_study/several_K_data_generation.R | 2 +- .../several_K_data_generation_uncertainties.R | 2 +- simulations_study/several_K_exploitation.R | 4 +- .../several_K_exploitation_FUN.R | 4 +- .../simulation_exploitation_bayou_design.R | 2 +- ...on_exploitation_bayou_design_alpha_known.R | 2 +- tests/testthat/test-core-log_likelihood.R | 4 +- 29 files changed, 138 insertions(+), 138 deletions(-) rename man/{extract.variance_covariance.Rd => extract_variance_covariance.Rd} (78%) rename man/{prod.index.Rd => prod_index.Rd} (77%) rename man/{sum.partitions.Rd => sum_partitions.Rd} (78%) rename man/{sum.partitions.cardFixed.Rd => sum_partitions.cardFixed.Rd} (78%) rename man/{sum.prod.comb.Rd => sum_prod.comb.Rd} (79%) rename man/{sum.simplex.Rd => sum_simplex.Rd} (78%) diff --git a/R/E_step.R b/R/E_step.R index 419da49..cb8375c 100644 --- a/R/E_step.R +++ b/R/E_step.R @@ -29,7 +29,7 @@ #' @details #' This function takes parameters sim, Sigma and Sigma_YY_inv from #' \code{compute_mean_variance.simple}. It uses functions -#' \code{extract.variance_covariance}, \code{extract.covariance_parents}, and +#' \code{extract_variance_covariance}, \code{extract_covariance_parents}, and #' \code{extract_simulate_internal} to extract the needed quantities from these objects. #' #' @param phylo Input tree. @@ -126,8 +126,8 @@ compute_cond_law.simple <- function (phylo, Y_data_vec, sim, where="nodes", what="optimal.values")) # NULL if BM ## Variance Covariance - Sigma_YZ <- extract.variance_covariance(Sigma, what="YZ", masque_data) - Sigma_ZZ <- extract.variance_covariance(Sigma, what="ZZ", masque_data) + Sigma_YZ <- extract_variance_covariance(Sigma, what="YZ", masque_data) + Sigma_ZZ <- extract_variance_covariance(Sigma, what="ZZ", masque_data) # temp <- Sigma_YZ %*% Sigma_YY_inv temp <- Sigma_YZ %*% Sigma_YY_chol_inv # Y_data_vec <- as.vector(Y_data) @@ -173,12 +173,12 @@ compute_cond_law.simple <- function (phylo, Y_data_vec, sim, } } # Nodes - varariances - var_nodes <- extract.variance_nodes(phylo, + var_nodes <- extract_variance_nodes(phylo, conditional_variance_covariance_nodes) conditional_law_X$variances <- array(c(var_tips, var_nodes), c(p, p, ntaxa + Nnode)) # Nodes - covariances - cov_nodes <- extract.covariance_parents(phylo, + cov_nodes <- extract_covariance_parents(phylo, conditional_variance_covariance_nodes) conditional_law_X$covariances <- array(c(cov_tips, cov_nodes), c(p, p, ntaxa + Nnode)) @@ -290,9 +290,9 @@ compute_cond_law.simple.nomissing.BM <- function (phylo, Y_data, sim, ## compute_fixed_moments <- function(times_shared, ntaxa){ masque_data <- c(rep(TRUE, ntaxa), rep(FALSE, dim(times_shared)[1] - ntaxa)) - C_YZ <- extract.variance_covariance(times_shared, what="YZ", masque_data) - C_YY <- extract.variance_covariance(times_shared, what="YY", masque_data) - C_ZZ <- extract.variance_covariance(times_shared, what="ZZ", masque_data) + C_YZ <- extract_variance_covariance(times_shared, what="YZ", masque_data) + C_YY <- extract_variance_covariance(times_shared, what="YY", masque_data) + C_ZZ <- extract_variance_covariance(times_shared, what="ZZ", masque_data) C_YY_chol <- chol(C_YY) C_YY_chol_inv <- backsolve(C_YY_chol, diag(ncol(C_YY_chol))) temp <- C_YZ %*% C_YY_chol_inv @@ -310,7 +310,7 @@ compute_fixed_moments <- function(times_shared, ntaxa){ #' @title Extract sub-matrices of variance. #' #' @description -#' \code{extract.variance_covariance} return the adequate sub-matrix. +#' \code{extract_variance_covariance} return the adequate sub-matrix. #' #' @param struct structural matrix of size (ntaxa+Nnode)*p, result #' of function \code{compute_variance_covariance} @@ -318,14 +318,14 @@ compute_fixed_moments <- function(times_shared, ntaxa){ #' "YY" : sub-matrix of tips (p*ntaxa first lines and columns) #' "YZ" : sub matrix tips x nodes (p*Nnode last rows and p*ntaxa first columns) #' "ZZ" : sub matrix of nodes (p*Nnode last rows and columns) -#' @param miss missing values of Y_data +#' @param masque_data Mask of missing data #' #' @return sub-matrix of variance covariance. #' #' @keywords internal #' ## -extract.variance_covariance <- function(struct, what=c("YY","YZ","ZZ"), +extract_variance_covariance <- function(struct, what=c("YY","YZ","ZZ"), masque_data = c(rep(TRUE, attr(struct, "ntaxa") * attr(struct, "p_dim")), rep(FALSE, (dim(struct)[1] - attr(struct, "ntaxa")) * attr(struct, "p_dim")))){ # ntaxa <- attr(struct, "ntaxa") @@ -343,7 +343,7 @@ extract.variance_covariance <- function(struct, what=c("YY","YZ","ZZ"), } ## -# extract.covariance_parents (phylo, struct) +# extract_covariance_parents (phylo, struct) # PARAMETERS: # @phylo (tree) # @struct (matrix) structural matrix of size ntaxa+Nnode, result of function compute_times_ca, compute_dist_phy or compute_variance_covariance @@ -358,7 +358,7 @@ extract.variance_covariance <- function(struct, what=c("YY","YZ","ZZ"), # REVISIONS: # 22/05/14 - Initial release ## -extract.covariance_parents <- function(phylo, struct){ +extract_covariance_parents <- function(phylo, struct){ ntaxa <- length(phylo$tip.label) p <- attr(struct, "p_dim") m <- dim(phylo$edge)[1] - ntaxa + 1 @@ -391,7 +391,7 @@ extract.covariance_parents <- function(phylo, struct){ return(arr) } -extract.variance_nodes <- function(phylo, struct){ +extract_variance_nodes <- function(phylo, struct){ ntaxa <- length(phylo$tip.label) p <- attr(struct, "p_dim") m <- dim(phylo$edge)[1] - ntaxa + 1 @@ -857,7 +857,7 @@ compute_mean_variance.simple <- function(phylo, Sigma <- compute_variance_covariance(times_shared = times_shared, distances_phylo = distances_phylo, params_old = params_old) - Sigma_YY <- extract.variance_covariance(Sigma, what="YY", masque_data = masque_data) + Sigma_YY <- extract_variance_covariance(Sigma, what="YY", masque_data = masque_data) Sigma_YY_chol <- chol(Sigma_YY) Sigma_YY_chol_inv <- backsolve(Sigma_YY_chol, Matrix::diag(ncol(Sigma_YY_chol))) #Sigma_YY_inv <- chol2inv(Sigma_YY_chol) @@ -975,7 +975,7 @@ compute_mahalanobis_distance.simple.nomissing.BM <- function(phylo, Y_data, sim, #' @details #' This function takes parameters sim, Sigma and Sigma_YY_inv from #' \code{compute_mean_variance.simple}. It uses functions -#' \code{extract.variance_covariance}, \code{extract.covariance_parents}, and +#' \code{extract_variance_covariance}, \code{extract_covariance_parents}, and #' \code{extract_simulate_internal} to extract the needed quantities from these objects. #' #' @param phylo Input tree. @@ -994,7 +994,7 @@ compute_log_likelihood.simple <- function(phylo, Y_data_vec, sim, masque_data = c(rep(TRUE, dim(sim)[1] * length(phylo$tip.label)), rep(FALSE, dim(sim)[1] * phylo$Nnode)), ...){ # ntaxa <- length(phylo$tip.label) - Sigma_YY <- extract.variance_covariance(Sigma, what="YY", masque_data) + Sigma_YY <- extract_variance_covariance(Sigma, what="YY", masque_data) logdetSigma_YY <- Matrix::determinant(Sigma_YY, logarithm = TRUE)$modulus m_Y <- extract_simulate_internal(sim, where="tips", what="expectations") LL <- length(Y_data_vec) * log(2*pi) + logdetSigma_YY @@ -1020,7 +1020,7 @@ compute_log_likelihood.simple.nomissing.BM <- function(phylo, Y_data, sim, # Sigma, Sigma_YY_chol_inv, # miss, masque_data){ # ntaxa <- length(phylo$tip.label) -# Sigma_YY <- extract.variance_covariance(Sigma, what="YY", masque_data) +# Sigma_YY <- extract_variance_covariance(Sigma, what="YY", masque_data) # logdetSigma_YY <- determinant(Sigma_YY, logarithm = TRUE)$modulus # return( - (ntaxa * log(2*pi) + logdetSigma_YY) / 2) # } @@ -1036,8 +1036,8 @@ compute_log_likelihood.simple.nomissing.BM <- function(phylo, Y_data, sim, # } # compute_entropy.simple <- function(Sigma, Sigma_YY_inv){ -# Sigma_YZ <- extract.variance_covariance(Sigma, what="YZ") -# Sigma_ZZ <- extract.variance_covariance(Sigma, what="ZZ") +# Sigma_YZ <- extract_variance_covariance(Sigma, what="YZ") +# Sigma_ZZ <- extract_variance_covariance(Sigma, what="ZZ") # conditional_variance_covariance <- Sigma_ZZ - Sigma_YZ%*%Sigma_YY_inv%*%t(Sigma_YZ) # N <- dim(Sigma_ZZ)[1] # logdet_conditional_variance_covariance <- determinant(conditional_variance_covariance, logarithm = TRUE)$modulus diff --git a/R/M_step.R b/R/M_step.R index 97c438a..c391f95 100644 --- a/R/M_step.R +++ b/R/M_step.R @@ -543,7 +543,7 @@ compute_var_diff.OU <- function(phylo, conditional_law_X, selection.strength) { #' @title Compute weighted sum of var_diff #' #' @description -#' \code{compute_sum_var_diff} computes \sum_\{e edge\} ell_j * Var[X_j - X_pa(j) | Y] +#' \code{compute_sum_var_diff} computes sum_\{e edge\} ell_j * Var[X_j - X_pa(j) | Y] #' #' @param phylo a phylogenetic tree #' @param var_diff result of function \code{compute_var_diff.BM} diff --git a/R/init_EM.R b/R/init_EM.R index 38549aa..9bbb9e2 100644 --- a/R/init_EM.R +++ b/R/init_EM.R @@ -1664,7 +1664,7 @@ init.EM.lasso <- function(phylo, Fm <- compute_tree_correlations_matrix(times_shared = times_shared, distances_phylo = distances_phylo, params_old = params_sigma) - Fm_YY <- extract.variance_covariance(Fm, what="YY", + Fm_YY <- extract_variance_covariance(Fm, what="YY", masque_data = c(rep(TRUE, ntaxa), rep(FALSE, dim(Fm)[1] - ntaxa))) # Cholesky of tree-correlations diff --git a/R/model_selection.R b/R/model_selection.R index 19fc4af..e7dde61 100644 --- a/R/model_selection.R +++ b/R/model_selection.R @@ -348,7 +348,7 @@ penalty_pBIC_scalarOU <- function(params, tree, times_shared, OU = compute_tree_correlations_matrix.scOU, scOU = compute_tree_correlations_matrix.scOU) C <- compute_tree_correlations_matrix(times_shared, distances_phylo, params) - C <- extract.variance_covariance(C, what="YY", + C <- extract_variance_covariance(C, what="YY", masque_data = c(rep(TRUE, ntaxa), rep(FALSE, dim(C)[1] - ntaxa))) C <- 1/(2*as.vector(params$selection.strength)) * C @@ -382,7 +382,7 @@ penalty_pBIC_independent <- function(params, model_complexity, tree, OU = compute_tree_correlations_matrix.scOU, scOU = compute_tree_correlations_matrix.scOU) C <- compute_tree_correlations_matrix(times_shared, distances_phylo, params) - C <- extract.variance_covariance(C, what="YY", + C <- extract_variance_covariance(C, what="YY", masque_data = c(rep(TRUE, ntaxa), rep(FALSE, dim(C)[1] - ntaxa))) C_inv <- solve(C) @@ -453,7 +453,7 @@ penalty_pBIC_l1ou_unit <- function(params, tree, times_shared, distances_phylo, OU = compute_tree_correlations_matrix.scOU, scOU = compute_tree_correlations_matrix.scOU) C <- compute_tree_correlations_matrix(times_shared, distances_phylo, params) - C <- extract.variance_covariance(C, what="YY", + C <- extract_variance_covariance(C, what="YY", masque_data = c(rep(TRUE, ntaxa), rep(FALSE, dim(C)[1] - ntaxa))) C <- 1/(2*as.vector(params$selection.strength)) * C diff --git a/R/partitionsNumber.R b/R/partitionsNumber.R index e60cd1f..12f43eb 100644 --- a/R/partitionsNumber.R +++ b/R/partitionsNumber.R @@ -222,7 +222,7 @@ update.partitionsNumber.bin <- function(daughtersParams, ...){ #' Nk and Ak of a node given its daughters. #' #' @details -#' Uses functions \code{sum.partitions} and \code{sum.simplex} to compute +#' Uses functions \code{sum_partitions} and \code{sum_simplex} to compute #' the needed sums. #' #' @param daughtersParams : matrix with 2*npart columns. Each row contains @@ -250,8 +250,8 @@ update.partitionsNumber.gen <- function(daughtersParams, ...){ for (K in 2:npart){ N <- N_daughters[,1:K] A <- A_daughters[,1:K] - nbrComp[K] <- sum.simplex(N, K, p) + sum.partitions(A, N, K, p, 2) - nbrMarkedComp[K] <- sum.partitions(A, N, K, p, 1) + nbrComp[K] <- sum_simplex(N, K, p) + sum_partitions(A, N, K, p, 2) + nbrMarkedComp[K] <- sum_partitions(A, N, K, p, 1) } nbrAdm <- c(nbrComp, nbrMarkedComp) return(nbrAdm) @@ -284,10 +284,10 @@ update.partitionsNumber.gen <- function(daughtersParams, ...){ #' @title Product of elements of a matrix #' #' @description -#' \code{prod.index} return the product of chosen elements of a matrix. +#' \code{prod_index} return the product of chosen elements of a matrix. #' #' @details -#' This function is to be used in \code{sum.simplex} to be applied to all the +#' This function is to be used in \code{sum_simplex} to be applied to all the #' elements given by \code{xsimplex}. #' Performs the product : X[1,Id[1]]*X[2,Id[2]]*...*X[p,Id[p]] if all the #' elements of Id are positive. Otherwise, just return 0. @@ -302,7 +302,7 @@ update.partitionsNumber.gen <- function(daughtersParams, ...){ #' #05/06/14 - Initial release ## -prod.index <- function(X,Id){ +prod_index <- function(X,Id){ if (0 %in% Id) return(0) # Indice hors limites return(prod(diag(X[,Id]))) # if (any(Id == 0)) return(0) # Indice hors limites @@ -313,13 +313,13 @@ prod.index <- function(X,Id){ #' @title Sum on a simplex #' #' @description -#' \code{sum.simplex} returns the sum on k_1+...+k_p=K, k_i>0 of the products +#' \code{sum_simplex} returns the sum on k_1+...+k_p=K, k_i>0 of the products #' of NN[i,k_i]. #' #' @details #' This function uses \code{xsimplex} to perform the product of NN[i,k_i] for #' all combination of k_i such that k_1+...+k_p=K, k_i>0, using function -#' \code{prod.index}. Then sum all the products. +#' \code{prod_index}. Then sum all the products. #' #' @param NN a matrix with p rows and K column. Each row contains the number of #' partition in 1<=k<=K groups for one of the p children of a given node @@ -332,8 +332,8 @@ prod.index <- function(X,Id){ #' #05/06/14 - Initial release ## -sum.simplex <- function (NN, K, p) { - return(sum(combinat::xsimplex(p, K, fun = prod.index, +sum_simplex <- function (NN, K, p) { + return(sum(combinat::xsimplex(p, K, fun = prod_index, simplify = TRUE, X = NN))) } @@ -342,11 +342,11 @@ sum.simplex <- function (NN, K, p) { #' @title Sum on a simplex #' #' @description -#' \code{sum.prod.comb} returns the sum on k_1+...+k_p=K+|I|-1, k_i>0 of the +#' \code{sum_prod.comb} returns the sum on k_1+...+k_p=K+|I|-1, k_i>0 of the #' products of prod(A[i,k_i], i in I)*prod(N[i,k_i], i not in I). #' #' @details -#' This function uses \code{sum.simplex} to perform the wanted sum on a ad-hoc +#' This function uses \code{sum_simplex} to perform the wanted sum on a ad-hoc #' matrix, combination of rows of A and N. #' #' @param I a vector of integers representing a subset of [1,p] @@ -364,20 +364,20 @@ sum.simplex <- function (NN, K, p) { #' #05/06/14 - Initial release ## -sum.prod.comb <- function(I, A, N, K, p){ +sum_prod.comb <- function(I, A, N, K, p){ KK <- K + length(I) - 1 AN <- matrix(NA, nrow = p, ncol = K) AN[I, ] <- A[I, ] AN[-I, ] <- N[-I, ] - return(sum.simplex(AN, KK, p)) + return(sum_simplex(AN, KK, p)) } ## #' @title Sum on subsets of a given cardinal. #' #' @description -#' \code{sum.partitions.cardFixed} returns the sum on I subset of [1,p], |I| -#' fixed, of the sums computed by \code{sum.prod.comb}. +#' \code{sum_partitions.cardFixed} returns the sum on I subset of [1,p], |I| +#' fixed, of the sums computed by \code{sum_prod.comb}. #' #' @details #' This function uses \code{combn} to enumerate all the subsets I of [1,p] of @@ -398,9 +398,9 @@ sum.prod.comb <- function(I, A, N, K, p){ #' #05/06/14 - Initial release ## -sum.partitions.cardFixed <- function(A, N, K, p, cardI){ +sum_partitions.cardFixed <- function(A, N, K, p, cardI){ return(sum(combinat::combn(p, cardI, - fun = sum.prod.comb, + fun = sum_prod.comb, simplify = TRUE, A = A, N = N, K = K, p = p))) } @@ -409,10 +409,10 @@ sum.partitions.cardFixed <- function(A, N, K, p, cardI){ #' @title Sum on all subsets. #' #' @description -#' \code{sum.partitions} returns the sum on all I subset of [1,p], with |I|>m. +#' \code{sum_partitions} returns the sum on all I subset of [1,p], with |I|>m. #' #' @details -#' This function applies \code{sum.partitions.cardFixed} to all integer between +#' This function applies \code{sum_partitions.cardFixed} to all integer between #' m and p, and sum the results. #' #' @param A a matrix with p rows and K column. Each row contains the number of @@ -430,8 +430,8 @@ sum.partitions.cardFixed <- function(A, N, K, p, cardI){ #' #05/06/14 - Initial release ## -sum.partitions <- function(A, N, K, p, m) { - return(sum(sapply(m:p, function(x) sum.partitions.cardFixed(A, N, K, p, x)))) +sum_partitions <- function(A, N, K, p, m) { + return(sum(sapply(m:p, function(x) sum_partitions.cardFixed(A, N, K, p, x)))) } ############################################################################### diff --git a/R/plot_functions.R b/R/plot_functions.R index f6588fe..a934979 100644 --- a/R/plot_functions.R +++ b/R/plot_functions.R @@ -348,7 +348,7 @@ plot.PhyloEM <- function(x, if (missing(imposed_scale)) imposed_scale <- Y_state ## Plotting - plot.data.process.actual(Y.state = Y_state, + plot_data.process.actual(Y.state = Y_state, phylo = x$phylo, params = params, process = x$process, @@ -390,7 +390,7 @@ plot.PhyloEM <- function(x, ...) } -plot.data.process.actual <- function(Y.state, phylo, params, +plot_data.process.actual <- function(Y.state, phylo, params, miss = is.na(Y.state), process = "BM", #norm = max(abs(Y.state)), diff --git a/R/simulate.R b/R/simulate.R index 47de8fa..039e9a1 100644 --- a/R/simulate.R +++ b/R/simulate.R @@ -245,7 +245,7 @@ plot.params_process <- function(x, if (missing(imposed_scale)) imposed_scale <- Y_state ## Plotting - plot.data.process.actual(Y.state = Y_state[traits, , drop = FALSE], + plot_data.process.actual(Y.state = Y_state[traits, , drop = FALSE], phylo = phylo, params = params, process = x$process, diff --git a/man/compute_E.simple.Rd b/man/compute_E.simple.Rd index ede7b62..c0fcc20 100644 --- a/man/compute_E.simple.Rd +++ b/man/compute_E.simple.Rd @@ -43,7 +43,7 @@ values beta(t_j) \details{ This function takes parameters sim, Sigma and Sigma_YY_inv from \code{compute_mean_variance.simple}. It uses functions -\code{extract.variance_covariance}, \code{extract.covariance_parents}, and +\code{extract_variance_covariance}, \code{extract_covariance_parents}, and \code{extract_simulate_internal} to extract the needed quantities from these objects. } \keyword{internal} diff --git a/man/compute_log_likelihood.simple.Rd b/man/compute_log_likelihood.simple.Rd index c336555..fc1c6a2 100644 --- a/man/compute_log_likelihood.simple.Rd +++ b/man/compute_log_likelihood.simple.Rd @@ -31,7 +31,7 @@ in the simple case where the inverse of the variance matrix is given. \details{ This function takes parameters sim, Sigma and Sigma_YY_inv from \code{compute_mean_variance.simple}. It uses functions -\code{extract.variance_covariance}, \code{extract.covariance_parents}, and +\code{extract_variance_covariance}, \code{extract_covariance_parents}, and \code{extract_simulate_internal} to extract the needed quantities from these objects. } \keyword{internal} diff --git a/man/compute_sum_var_diff.Rd b/man/compute_sum_var_diff.Rd index e2bd43a..e49e44f 100644 --- a/man/compute_sum_var_diff.Rd +++ b/man/compute_sum_var_diff.Rd @@ -15,6 +15,6 @@ compute_sum_var_diff(phylo, var_diff) matrix p x p } \description{ -\code{compute_sum_var_diff} computes \sum_\{e edge\} ell_j * Var[X_j - X_pa(j) | Y] +\code{compute_sum_var_diff} computes sum_\{e edge\} ell_j * Var[X_j - X_pa(j) | Y] } \keyword{internal} diff --git a/man/extract.variance_covariance.Rd b/man/extract_variance_covariance.Rd similarity index 78% rename from man/extract.variance_covariance.Rd rename to man/extract_variance_covariance.Rd index 88e84f2..690dc43 100644 --- a/man/extract.variance_covariance.Rd +++ b/man/extract_variance_covariance.Rd @@ -1,10 +1,10 @@ % Generated by roxygen2: do not edit by hand % Please edit documentation in R/E_step.R -\name{extract.variance_covariance} -\alias{extract.variance_covariance} +\name{extract_variance_covariance} +\alias{extract_variance_covariance} \title{Extract sub-matrices of variance.} \usage{ -\method{extract}{variance_covariance}( +extract_variance_covariance( struct, what = c("YY", "YZ", "ZZ"), masque_data = c(rep(TRUE, attr(struct, "ntaxa") * attr(struct, "p_dim")), rep(FALSE, @@ -20,12 +20,12 @@ of function \code{compute_variance_covariance}} "YZ" : sub matrix tips x nodes (p*Nnode last rows and p*ntaxa first columns) "ZZ" : sub matrix of nodes (p*Nnode last rows and columns)} -\item{miss}{missing values of Y_data} +\item{masque_data}{Mask of missing data} } \value{ sub-matrix of variance covariance. } \description{ -\code{extract.variance_covariance} return the adequate sub-matrix. +\code{extract_variance_covariance} return the adequate sub-matrix. } \keyword{internal} diff --git a/man/prod.index.Rd b/man/prod_index.Rd similarity index 77% rename from man/prod.index.Rd rename to man/prod_index.Rd index 07bb084..aa26d6c 100644 --- a/man/prod.index.Rd +++ b/man/prod_index.Rd @@ -1,10 +1,10 @@ % Generated by roxygen2: do not edit by hand % Please edit documentation in R/partitionsNumber.R -\name{prod.index} -\alias{prod.index} +\name{prod_index} +\alias{prod_index} \title{Product of elements of a matrix} \usage{ -\method{prod}{index}(X, Id) +prod_index(X, Id) } \arguments{ \item{X}{a matrix with p rows and K column. Each row contains the number of @@ -16,10 +16,10 @@ partition in 1<=k<=K groups for one of the p children of a given node} double : the result of the product. } \description{ -\code{prod.index} return the product of chosen elements of a matrix. +\code{prod_index} return the product of chosen elements of a matrix. } \details{ -This function is to be used in \code{sum.simplex} to be applied to all the +This function is to be used in \code{sum_simplex} to be applied to all the elements given by \code{xsimplex}. Performs the product : X[1,Id[1]]*X[2,Id[2]]*...*X[p,Id[p]] if all the elements of Id are positive. Otherwise, just return 0. diff --git a/man/sum.partitions.Rd b/man/sum_partitions.Rd similarity index 78% rename from man/sum.partitions.Rd rename to man/sum_partitions.Rd index a840444..f2ac0d8 100644 --- a/man/sum.partitions.Rd +++ b/man/sum_partitions.Rd @@ -1,10 +1,10 @@ % Generated by roxygen2: do not edit by hand % Please edit documentation in R/partitionsNumber.R -\name{sum.partitions} -\alias{sum.partitions} +\name{sum_partitions} +\alias{sum_partitions} \title{Sum on all subsets.} \usage{ -\method{sum}{partitions}(A, N, K, p, m) +sum_partitions(A, N, K, p, m) } \arguments{ \item{A}{a matrix with p rows and K column. Each row contains the number of @@ -24,10 +24,10 @@ partition in 1<=k<=K groups for one of the p children of a given node} double : the result of the sum. } \description{ -\code{sum.partitions} returns the sum on all I subset of [1,p], with |I|>m. +\code{sum_partitions} returns the sum on all I subset of [1,p], with |I|>m. } \details{ -This function applies \code{sum.partitions.cardFixed} to all integer between +This function applies \code{sum_partitions.cardFixed} to all integer between m and p, and sum the results. } \keyword{internal} diff --git a/man/sum.partitions.cardFixed.Rd b/man/sum_partitions.cardFixed.Rd similarity index 78% rename from man/sum.partitions.cardFixed.Rd rename to man/sum_partitions.cardFixed.Rd index 218f85b..778f556 100644 --- a/man/sum.partitions.cardFixed.Rd +++ b/man/sum_partitions.cardFixed.Rd @@ -1,10 +1,10 @@ % Generated by roxygen2: do not edit by hand % Please edit documentation in R/partitionsNumber.R -\name{sum.partitions.cardFixed} -\alias{sum.partitions.cardFixed} +\name{sum_partitions.cardFixed} +\alias{sum_partitions.cardFixed} \title{Sum on subsets of a given cardinal.} \usage{ -\method{sum}{partitions.cardFixed}(A, N, K, p, cardI) +sum_partitions.cardFixed(A, N, K, p, cardI) } \arguments{ \item{A}{a matrix with p rows and K column. Each row contains the number of @@ -24,8 +24,8 @@ partition in 1<=k<=K groups for one of the p children of a given node} double : the result of the sum. } \description{ -\code{sum.partitions.cardFixed} returns the sum on I subset of [1,p], |I| -fixed, of the sums computed by \code{sum.prod.comb}. +\code{sum_partitions.cardFixed} returns the sum on I subset of [1,p], |I| +fixed, of the sums computed by \code{sum_prod.comb}. } \details{ This function uses \code{combn} to enumerate all the subsets I of [1,p] of diff --git a/man/sum.prod.comb.Rd b/man/sum_prod.comb.Rd similarity index 79% rename from man/sum.prod.comb.Rd rename to man/sum_prod.comb.Rd index d695244..17f8121 100644 --- a/man/sum.prod.comb.Rd +++ b/man/sum_prod.comb.Rd @@ -1,10 +1,10 @@ % Generated by roxygen2: do not edit by hand % Please edit documentation in R/partitionsNumber.R -\name{sum.prod.comb} -\alias{sum.prod.comb} +\name{sum_prod.comb} +\alias{sum_prod.comb} \title{Sum on a simplex} \usage{ -\method{sum}{prod.comb}(I, A, N, K, p) +sum_prod.comb(I, A, N, K, p) } \arguments{ \item{I}{a vector of integers representing a subset of [1,p]} @@ -24,11 +24,11 @@ partition in 1<=k<=K groups for one of the p children of a given node} double : the result of the sum. } \description{ -\code{sum.prod.comb} returns the sum on k_1+...+k_p=K+|I|-1, k_i>0 of the +\code{sum_prod.comb} returns the sum on k_1+...+k_p=K+|I|-1, k_i>0 of the products of prod(A[i,k_i], i in I)*prod(N[i,k_i], i not in I). } \details{ -This function uses \code{sum.simplex} to perform the wanted sum on a ad-hoc +This function uses \code{sum_simplex} to perform the wanted sum on a ad-hoc matrix, combination of rows of A and N. } \keyword{internal} diff --git a/man/sum.simplex.Rd b/man/sum_simplex.Rd similarity index 78% rename from man/sum.simplex.Rd rename to man/sum_simplex.Rd index 8987a73..5b1fce8 100644 --- a/man/sum.simplex.Rd +++ b/man/sum_simplex.Rd @@ -1,10 +1,10 @@ % Generated by roxygen2: do not edit by hand % Please edit documentation in R/partitionsNumber.R -\name{sum.simplex} -\alias{sum.simplex} +\name{sum_simplex} +\alias{sum_simplex} \title{Sum on a simplex} \usage{ -\method{sum}{simplex}(NN, K, p) +sum_simplex(NN, K, p) } \arguments{ \item{NN}{a matrix with p rows and K column. Each row contains the number of @@ -18,12 +18,12 @@ partition in 1<=k<=K groups for one of the p children of a given node} double : the result of the sum. } \description{ -\code{sum.simplex} returns the sum on k_1+...+k_p=K, k_i>0 of the products +\code{sum_simplex} returns the sum on k_1+...+k_p=K, k_i>0 of the products of NN[i,k_i]. } \details{ This function uses \code{xsimplex} to perform the product of NN[i,k_i] for all combination of k_i such that k_1+...+k_p=K, k_i>0, using function -\code{prod.index}. Then sum all the products. +\code{prod_index}. Then sum all the products. } \keyword{internal} diff --git a/man/update.partitionsNumber.gen.Rd b/man/update.partitionsNumber.gen.Rd index 43ce8fb..4c594f3 100644 --- a/man/update.partitionsNumber.gen.Rd +++ b/man/update.partitionsNumber.gen.Rd @@ -22,7 +22,7 @@ starting at the current node. Nk and Ak of a node given its daughters. } \details{ -Uses functions \code{sum.partitions} and \code{sum.simplex} to compute +Uses functions \code{sum_partitions} and \code{sum_simplex} to compute the needed sums. } \keyword{internal} diff --git a/simulations_study/multivariate_data_generation.R b/simulations_study/multivariate_data_generation.R index d5dc0d1..6cbd509 100644 --- a/simulations_study/multivariate_data_generation.R +++ b/simulations_study/multivariate_data_generation.R @@ -169,7 +169,7 @@ W <- compute_actualization_matrix_ultrametric(trees[["128"]], alpha_base * diag( vec_Y <- kronecker(T_tree[["128"]], diag(1, p_base, p_base)) %*% W %*% as.vector(Delta) X1.tips.exp.mat <- matrix(vec_Y, p_base, 128) + beta_0 unique(X1.tips.exp.mat[1, ]) -plot.data.process.actual(Y.state = X1.tips.exp.mat[1, ], +plot_data.process.actual(Y.state = X1.tips.exp.mat[1, ], phylo = trees[["128"]], params = list(shifts = list(edges = shifts_grid[["128_3"]]$edges, values = shifts_grid[["128_3"]]$values[1, ])), adj.root = 0, @@ -199,7 +199,7 @@ extract.parsimonyNumber(parsimonyNumber(trees[["128"]], clusters_from_shifts(tre # vec_Y <- kronecker(T_tree[["128"]], diag(1, p_base, p_base)) %*% W %*% as.vector(Delta) # X1.tips.exp.mat <- matrix(vec_Y, p_base, ntaxa_base) + beta_0 # unique(X1.tips.exp.mat[1, ]) -# plot.data.process.actual(Y.state = X1.tips.exp.mat[1, ], +# plot_data.process.actual(Y.state = X1.tips.exp.mat[1, ], # phylo = trees[["128"]], # params = list(shifts = list(edges = shifts_grid[["128_7"]]$edges, values = shifts_grid[["128_7"]]$values[1, ])), # adj.root = 0, @@ -234,7 +234,7 @@ extract.parsimonyNumber(parsimonyNumber(trees[["128"]], clusters_from_shifts(tre # vec_Y <- kronecker(T_tree[["128"]], diag(1, p_base, p_base)) %*% W %*% as.vector(Delta) # X1.tips.exp.mat <- matrix(vec_Y, p_base, ntaxa_base) + beta_0 # unique(X1.tips.exp.mat[1, ]) -# plot.data.process.actual(Y.state = X1.tips.exp.mat[1, ], +# plot_data.process.actual(Y.state = X1.tips.exp.mat[1, ], # phylo = trees[["128"]], # params = list(shifts = list(edges = shifts_grid[["128_11"]]$edges, values = shifts_grid[["128_11"]]$values[1, ])), # adj.root = 0, @@ -274,7 +274,7 @@ extract.parsimonyNumber(parsimonyNumber(trees[["128"]], clusters_from_shifts(tre # vec_Y <- kronecker(T_tree[["128"]], diag(1, p_base, p_base)) %*% W %*% as.vector(Delta) # X1.tips.exp.mat <- matrix(vec_Y, p_base, ntaxa_base) + beta_0 # unique(X1.tips.exp.mat[1, ]) -# plot.data.process.actual(Y.state = X1.tips.exp.mat[1, ], +# plot_data.process.actual(Y.state = X1.tips.exp.mat[1, ], # phylo = trees[["128"]], # params = list(shifts = list(edges = shifts_grid[["128_15"]]$edges, values = shifts_grid[["128_15"]]$values[1, ])), # adj.root = 0, @@ -300,7 +300,7 @@ W <- compute_actualization_matrix_ultrametric(trees[["32"]], alpha_base * diag(1 vec_Y <- kronecker(T_tree[["32"]], diag(1, p_base, p_base)) %*% W %*% as.vector(Delta) X1.tips.exp.mat <- matrix(vec_Y, p_base, 32) + beta_0 unique(X1.tips.exp.mat[1, ]) -plot.data.process.actual(Y.state = X1.tips.exp.mat[1, ], +plot_data.process.actual(Y.state = X1.tips.exp.mat[1, ], phylo = trees[["32"]], params = list(shifts = list(edges = shifts_grid[["32_3"]]$edges, values = shifts_grid[["32_3"]]$values[1, ])), adj.root = 0, @@ -326,7 +326,7 @@ W <- compute_actualization_matrix_ultrametric(trees[["64"]], alpha_base * diag(1 vec_Y <- kronecker(T_tree[["64"]], diag(1, p_base, p_base)) %*% W %*% as.vector(Delta) X1.tips.exp.mat <- matrix(vec_Y, p_base, 64) + beta_0 unique(X1.tips.exp.mat[1, ]) -plot.data.process.actual(Y.state = X1.tips.exp.mat[1, ], +plot_data.process.actual(Y.state = X1.tips.exp.mat[1, ], phylo = trees[["64"]], params = list(shifts = list(edges = shifts_grid[["64_3"]]$edges, values = shifts_grid[["64_3"]]$values[1, ])), adj.root = 0, @@ -352,7 +352,7 @@ W <- compute_actualization_matrix_ultrametric(trees[["96"]], alpha_base * diag(1 vec_Y <- kronecker(T_tree[["96"]], diag(1, p_base, p_base)) %*% W %*% as.vector(Delta) X1.tips.exp.mat <- matrix(vec_Y, p_base, 96) + beta_0 unique(X1.tips.exp.mat[1, ]) -plot.data.process.actual(Y.state = X1.tips.exp.mat[1, ], +plot_data.process.actual(Y.state = X1.tips.exp.mat[1, ], phylo = trees[["96"]], params = list(shifts = list(edges = shifts_grid[["96_3"]]$edges, values = shifts_grid[["96_3"]]$values[1, ])), adj.root = 0, @@ -378,7 +378,7 @@ W <- compute_actualization_matrix_ultrametric(trees[["160"]], alpha_base * diag( vec_Y <- kronecker(T_tree[["160"]], diag(1, p_base, p_base)) %*% W %*% as.vector(Delta) X1.tips.exp.mat <- matrix(vec_Y, p_base, 160) + beta_0 unique(X1.tips.exp.mat[1, ]) -plot.data.process.actual(Y.state = X1.tips.exp.mat[1, ], +plot_data.process.actual(Y.state = X1.tips.exp.mat[1, ], phylo = trees[["160"]], params = list(shifts = list(edges = shifts_grid[["160_3"]]$edges, values = shifts_grid[["160_3"]]$values[1, ])), adj.root = 0, @@ -408,7 +408,7 @@ W <- compute_actualization_matrix_ultrametric(trees[["160"]], alpha_base * diag( vec_Y <- kronecker(T_tree[["160"]], diag(1, p_base, p_base)) %*% W %*% as.vector(Delta) X1.tips.exp.mat <- matrix(vec_Y, p_base, 160) + beta_0 unique(X1.tips.exp.mat[1, ]) -plot.data.process.actual(Y.state = X1.tips.exp.mat[1, ], +plot_data.process.actual(Y.state = X1.tips.exp.mat[1, ], phylo = trees[["160"]], params = list(shifts = list(edges = shifts_grid[["160_7"]]$edges, values = shifts_grid[["160_7"]]$values[1, ])), adj.root = 0, @@ -443,7 +443,7 @@ W <- compute_actualization_matrix_ultrametric(trees[["160"]], alpha_base * diag( vec_Y <- kronecker(T_tree[["160"]], diag(1, p_base, p_base)) %*% W %*% as.vector(Delta) X1.tips.exp.mat <- matrix(vec_Y, p_base, 160) + beta_0 unique(X1.tips.exp.mat[1, ]) -plot.data.process.actual(Y.state = X1.tips.exp.mat[1, ], +plot_data.process.actual(Y.state = X1.tips.exp.mat[1, ], phylo = trees[["160"]], params = list(shifts = list(edges = shifts_grid[["160_11"]]$edges, values = shifts_grid[["160_11"]]$values[1, ])), adj.root = 0, @@ -483,7 +483,7 @@ W <- compute_actualization_matrix_ultrametric(trees[["160"]], alpha_base * diag( vec_Y <- kronecker(T_tree[["160"]], diag(1, p_base, p_base)) %*% W %*% as.vector(Delta) X1.tips.exp.mat <- matrix(vec_Y, p_base, 160) + beta_0 unique(X1.tips.exp.mat[1, ]) -plot.data.process.actual(Y.state = X1.tips.exp.mat[1, ], +plot_data.process.actual(Y.state = X1.tips.exp.mat[1, ], phylo = trees[["160"]], params = list(shifts = list(edges = shifts_grid[["160_15"]]$edges, values = shifts_grid[["160_15"]]$values[1, ])), adj.root = 0, @@ -509,7 +509,7 @@ W <- compute_actualization_matrix_ultrametric(trees[["192"]], alpha_base * diag( vec_Y <- kronecker(T_tree[["192"]], diag(1, p_base, p_base)) %*% W %*% as.vector(Delta) X1.tips.exp.mat <- matrix(vec_Y, p_base, 192) + beta_0 unique(X1.tips.exp.mat[1, ]) -plot.data.process.actual(Y.state = X1.tips.exp.mat[1, ], +plot_data.process.actual(Y.state = X1.tips.exp.mat[1, ], phylo = trees[["192"]], params = list(shifts = list(edges = shifts_grid[["192_3"]]$edges, values = shifts_grid[["192_3"]]$values[1, ])), adj.root = 0, @@ -535,7 +535,7 @@ W <- compute_actualization_matrix_ultrametric(trees[["215"]], alpha_base * diag( vec_Y <- kronecker(T_tree[["215"]], diag(1, p_base, p_base)) %*% W %*% as.vector(Delta) X1.tips.exp.mat <- matrix(vec_Y, p_base, 215) + beta_0 unique(X1.tips.exp.mat[1, ]) -plot.data.process.actual(Y.state = X1.tips.exp.mat[1, ], +plot_data.process.actual(Y.state = X1.tips.exp.mat[1, ], phylo = trees[["215"]], params = list(shifts = list(edges = shifts_grid[["215_3"]]$edges, values = shifts_grid[["215_3"]]$values[1, ])), adj.root = 0, diff --git a/simulations_study/multivariate_data_generation_two_bases.R b/simulations_study/multivariate_data_generation_two_bases.R index bd98d2f..e9b6013 100644 --- a/simulations_study/multivariate_data_generation_two_bases.R +++ b/simulations_study/multivariate_data_generation_two_bases.R @@ -173,7 +173,7 @@ W <- compute_actualization_matrix_ultrametric(trees[["128"]], alpha_base * diag( vec_Y <- kronecker(T_tree[["128"]], diag(1, p_base, p_base)) %*% W %*% as.vector(Delta) X1.tips.exp.mat <- matrix(vec_Y, p_base, 128) + beta_0 unique(X1.tips.exp.mat[1, ]) -plot.data.process.actual(Y.state = X1.tips.exp.mat[1, ], +plot_data.process.actual(Y.state = X1.tips.exp.mat[1, ], phylo = trees[["128"]], params = list(shifts = list(edges = shifts_grid[["128_3"]]$edges, values = shifts_grid[["128_3"]]$values[1, ])), adj.root = 0, @@ -203,7 +203,7 @@ extract.parsimonyNumber(parsimonyNumber(trees[["128"]], clusters_from_shifts(tre # vec_Y <- kronecker(T_tree[["128"]], diag(1, p_base, p_base)) %*% W %*% as.vector(Delta) # X1.tips.exp.mat <- matrix(vec_Y, p_base, ntaxa_base) + beta_0 # unique(X1.tips.exp.mat[1, ]) -# plot.data.process.actual(Y.state = X1.tips.exp.mat[1, ], +# plot_data.process.actual(Y.state = X1.tips.exp.mat[1, ], # phylo = trees[["128"]], # params = list(shifts = list(edges = shifts_grid[["128_7"]]$edges, values = shifts_grid[["128_7"]]$values[1, ])), # adj.root = 0, @@ -238,7 +238,7 @@ extract.parsimonyNumber(parsimonyNumber(trees[["128"]], clusters_from_shifts(tre # vec_Y <- kronecker(T_tree[["128"]], diag(1, p_base, p_base)) %*% W %*% as.vector(Delta) # X1.tips.exp.mat <- matrix(vec_Y, p_base, ntaxa_base) + beta_0 # unique(X1.tips.exp.mat[1, ]) -# plot.data.process.actual(Y.state = X1.tips.exp.mat[1, ], +# plot_data.process.actual(Y.state = X1.tips.exp.mat[1, ], # phylo = trees[["128"]], # params = list(shifts = list(edges = shifts_grid[["128_11"]]$edges, values = shifts_grid[["128_11"]]$values[1, ])), # adj.root = 0, @@ -278,7 +278,7 @@ extract.parsimonyNumber(parsimonyNumber(trees[["128"]], clusters_from_shifts(tre # vec_Y <- kronecker(T_tree[["128"]], diag(1, p_base, p_base)) %*% W %*% as.vector(Delta) # X1.tips.exp.mat <- matrix(vec_Y, p_base, ntaxa_base) + beta_0 # unique(X1.tips.exp.mat[1, ]) -# plot.data.process.actual(Y.state = X1.tips.exp.mat[1, ], +# plot_data.process.actual(Y.state = X1.tips.exp.mat[1, ], # phylo = trees[["128"]], # params = list(shifts = list(edges = shifts_grid[["128_15"]]$edges, values = shifts_grid[["128_15"]]$values[1, ])), # adj.root = 0, @@ -304,7 +304,7 @@ W <- compute_actualization_matrix_ultrametric(trees[["32"]], alpha_base * diag(1 vec_Y <- kronecker(T_tree[["32"]], diag(1, p_base, p_base)) %*% W %*% as.vector(Delta) X1.tips.exp.mat <- matrix(vec_Y, p_base, 32) + beta_0 unique(X1.tips.exp.mat[1, ]) -plot.data.process.actual(Y.state = X1.tips.exp.mat[1, ], +plot_data.process.actual(Y.state = X1.tips.exp.mat[1, ], phylo = trees[["32"]], params = list(shifts = list(edges = shifts_grid[["32_3"]]$edges, values = shifts_grid[["32_3"]]$values[1, ])), adj.root = 0, @@ -330,7 +330,7 @@ W <- compute_actualization_matrix_ultrametric(trees[["64"]], alpha_base * diag(1 vec_Y <- kronecker(T_tree[["64"]], diag(1, p_base, p_base)) %*% W %*% as.vector(Delta) X1.tips.exp.mat <- matrix(vec_Y, p_base, 64) + beta_0 unique(X1.tips.exp.mat[1, ]) -plot.data.process.actual(Y.state = X1.tips.exp.mat[1, ], +plot_data.process.actual(Y.state = X1.tips.exp.mat[1, ], phylo = trees[["64"]], params = list(shifts = list(edges = shifts_grid[["64_3"]]$edges, values = shifts_grid[["64_3"]]$values[1, ])), adj.root = 0, @@ -356,7 +356,7 @@ W <- compute_actualization_matrix_ultrametric(trees[["96"]], alpha_base * diag(1 vec_Y <- kronecker(T_tree[["96"]], diag(1, p_base, p_base)) %*% W %*% as.vector(Delta) X1.tips.exp.mat <- matrix(vec_Y, p_base, 96) + beta_0 unique(X1.tips.exp.mat[1, ]) -plot.data.process.actual(Y.state = X1.tips.exp.mat[1, ], +plot_data.process.actual(Y.state = X1.tips.exp.mat[1, ], phylo = trees[["96"]], params = list(shifts = list(edges = shifts_grid[["96_3"]]$edges, values = shifts_grid[["96_3"]]$values[1, ])), adj.root = 0, @@ -382,7 +382,7 @@ W <- compute_actualization_matrix_ultrametric(trees[["160"]], alpha_base * diag( vec_Y <- kronecker(T_tree[["160"]], diag(1, p_base, p_base)) %*% W %*% as.vector(Delta) X1.tips.exp.mat <- matrix(vec_Y, p_base, 160) + beta_0 unique(X1.tips.exp.mat[1, ]) -plot.data.process.actual(Y.state = X1.tips.exp.mat[1, ], +plot_data.process.actual(Y.state = X1.tips.exp.mat[1, ], phylo = trees[["160"]], params = list(shifts = list(edges = shifts_grid[["160_3"]]$edges, values = shifts_grid[["160_3"]]$values[1, ])), adj.root = 0, @@ -412,7 +412,7 @@ W <- compute_actualization_matrix_ultrametric(trees[["160"]], alpha_base * diag( vec_Y <- kronecker(T_tree[["160"]], diag(1, p_base, p_base)) %*% W %*% as.vector(Delta) X1.tips.exp.mat <- matrix(vec_Y, p_base, 160) + beta_0 unique(X1.tips.exp.mat[1, ]) -plot.data.process.actual(Y.state = X1.tips.exp.mat[1, ], +plot_data.process.actual(Y.state = X1.tips.exp.mat[1, ], phylo = trees[["160"]], params = list(shifts = list(edges = shifts_grid[["160_7"]]$edges, values = shifts_grid[["160_7"]]$values[1, ])), adj.root = 0, @@ -447,7 +447,7 @@ W <- compute_actualization_matrix_ultrametric(trees[["160"]], alpha_base * diag( vec_Y <- kronecker(T_tree[["160"]], diag(1, p_base, p_base)) %*% W %*% as.vector(Delta) X1.tips.exp.mat <- matrix(vec_Y, p_base, 160) + beta_0 unique(X1.tips.exp.mat[1, ]) -plot.data.process.actual(Y.state = X1.tips.exp.mat[1, ], +plot_data.process.actual(Y.state = X1.tips.exp.mat[1, ], phylo = trees[["160"]], params = list(shifts = list(edges = shifts_grid[["160_11"]]$edges, values = shifts_grid[["160_11"]]$values[1, ])), adj.root = 0, @@ -487,7 +487,7 @@ W <- compute_actualization_matrix_ultrametric(trees[["160"]], alpha_base * diag( vec_Y <- kronecker(T_tree[["160"]], diag(1, p_base, p_base)) %*% W %*% as.vector(Delta) X1.tips.exp.mat <- matrix(vec_Y, p_base, 160) + beta_0 unique(X1.tips.exp.mat[1, ]) -plot.data.process.actual(Y.state = X1.tips.exp.mat[1, ], +plot_data.process.actual(Y.state = X1.tips.exp.mat[1, ], phylo = trees[["160"]], params = list(shifts = list(edges = shifts_grid[["160_15"]]$edges, values = shifts_grid[["160_15"]]$values[1, ])), adj.root = 0, @@ -513,7 +513,7 @@ W <- compute_actualization_matrix_ultrametric(trees[["192"]], alpha_base * diag( vec_Y <- kronecker(T_tree[["192"]], diag(1, p_base, p_base)) %*% W %*% as.vector(Delta) X1.tips.exp.mat <- matrix(vec_Y, p_base, 192) + beta_0 unique(X1.tips.exp.mat[1, ]) -plot.data.process.actual(Y.state = X1.tips.exp.mat[1, ], +plot_data.process.actual(Y.state = X1.tips.exp.mat[1, ], phylo = trees[["192"]], params = list(shifts = list(edges = shifts_grid[["192_3"]]$edges, values = shifts_grid[["192_3"]]$values[1, ])), adj.root = 0, @@ -539,7 +539,7 @@ W <- compute_actualization_matrix_ultrametric(trees[["256"]], alpha_base * diag( vec_Y <- kronecker(T_tree[["256"]], diag(1, p_base, p_base)) %*% W %*% as.vector(Delta) X1.tips.exp.mat <- matrix(vec_Y, p_base, 256) + beta_0 unique(X1.tips.exp.mat[1, ]) -plot.data.process.actual(Y.state = X1.tips.exp.mat[1, ], +plot_data.process.actual(Y.state = X1.tips.exp.mat[1, ], phylo = trees[["256"]], params = list(shifts = list(edges = shifts_grid[["256_3"]]$edges, values = shifts_grid[["256_3"]]$values[1, ])), adj.root = 0, diff --git a/simulations_study/multivariate_data_generation_two_bases_with errors.R b/simulations_study/multivariate_data_generation_two_bases_with errors.R index b6690c6..963e61a 100644 --- a/simulations_study/multivariate_data_generation_two_bases_with errors.R +++ b/simulations_study/multivariate_data_generation_two_bases_with errors.R @@ -209,7 +209,7 @@ W <- PhylogeneticEM:::compute_actualization_matrix_ultrametric(trees[["128"]], a vec_Y <- kronecker(T_tree[["128"]], diag(1, p_base, p_base)) %*% W %*% as.vector(Delta) X1.tips.exp.mat <- matrix(vec_Y, p_base, 128) + beta_0 unique(X1.tips.exp.mat[1, ]) -# plot.data.process.actual(Y.state = X1.tips.exp.mat[1, ], +# plot_data.process.actual(Y.state = X1.tips.exp.mat[1, ], # phylo = trees[["128"]], # params = list(shifts = list(edges = shifts_grid[["128_3"]]$edges, values = shifts_grid[["128_3"]]$values[1, ])), # adj.root = 0, @@ -239,7 +239,7 @@ PhylogeneticEM:::extract.parsimonyNumber(parsimonyNumber(trees[["128"]], cluster # vec_Y <- kronecker(T_tree[["128"]], diag(1, p_base, p_base)) %*% W %*% as.vector(Delta) # X1.tips.exp.mat <- matrix(vec_Y, p_base, ntaxa_base) + beta_0 # unique(X1.tips.exp.mat[1, ]) -# plot.data.process.actual(Y.state = X1.tips.exp.mat[1, ], +# plot_data.process.actual(Y.state = X1.tips.exp.mat[1, ], # phylo = trees[["128"]], # params = list(shifts = list(edges = shifts_grid[["128_7"]]$edges, values = shifts_grid[["128_7"]]$values[1, ])), # adj.root = 0, @@ -274,7 +274,7 @@ PhylogeneticEM:::extract.parsimonyNumber(parsimonyNumber(trees[["128"]], cluster # vec_Y <- kronecker(T_tree[["128"]], diag(1, p_base, p_base)) %*% W %*% as.vector(Delta) # X1.tips.exp.mat <- matrix(vec_Y, p_base, ntaxa_base) + beta_0 # unique(X1.tips.exp.mat[1, ]) -# plot.data.process.actual(Y.state = X1.tips.exp.mat[1, ], +# plot_data.process.actual(Y.state = X1.tips.exp.mat[1, ], # phylo = trees[["128"]], # params = list(shifts = list(edges = shifts_grid[["128_11"]]$edges, values = shifts_grid[["128_11"]]$values[1, ])), # adj.root = 0, @@ -314,7 +314,7 @@ PhylogeneticEM:::extract.parsimonyNumber(parsimonyNumber(trees[["128"]], cluster # vec_Y <- kronecker(T_tree[["128"]], diag(1, p_base, p_base)) %*% W %*% as.vector(Delta) # X1.tips.exp.mat <- matrix(vec_Y, p_base, ntaxa_base) + beta_0 # unique(X1.tips.exp.mat[1, ]) -# plot.data.process.actual(Y.state = X1.tips.exp.mat[1, ], +# plot_data.process.actual(Y.state = X1.tips.exp.mat[1, ], # phylo = trees[["128"]], # params = list(shifts = list(edges = shifts_grid[["128_15"]]$edges, values = shifts_grid[["128_15"]]$values[1, ])), # adj.root = 0, @@ -340,7 +340,7 @@ W <- PhylogeneticEM:::compute_actualization_matrix_ultrametric(trees[["32"]], al vec_Y <- kronecker(T_tree[["32"]], diag(1, p_base, p_base)) %*% W %*% as.vector(Delta) X1.tips.exp.mat <- matrix(vec_Y, p_base, 32) + beta_0 unique(X1.tips.exp.mat[1, ]) -# plot.data.process.actual(Y.state = X1.tips.exp.mat[1, ], +# plot_data.process.actual(Y.state = X1.tips.exp.mat[1, ], # phylo = trees[["32"]], # params = list(shifts = list(edges = shifts_grid[["32_3"]]$edges, values = shifts_grid[["32_3"]]$values[1, ])), # adj.root = 0, @@ -366,7 +366,7 @@ W <- PhylogeneticEM:::compute_actualization_matrix_ultrametric(trees[["64"]], al vec_Y <- kronecker(T_tree[["64"]], diag(1, p_base, p_base)) %*% W %*% as.vector(Delta) X1.tips.exp.mat <- matrix(vec_Y, p_base, 64) + beta_0 unique(X1.tips.exp.mat[1, ]) -# plot.data.process.actual(Y.state = X1.tips.exp.mat[1, ], +# plot_data.process.actual(Y.state = X1.tips.exp.mat[1, ], # phylo = trees[["64"]], # params = list(shifts = list(edges = shifts_grid[["64_3"]]$edges, values = shifts_grid[["64_3"]]$values[1, ])), # adj.root = 0, @@ -392,7 +392,7 @@ W <- PhylogeneticEM:::compute_actualization_matrix_ultrametric(trees[["96"]], al vec_Y <- kronecker(T_tree[["96"]], diag(1, p_base, p_base)) %*% W %*% as.vector(Delta) X1.tips.exp.mat <- matrix(vec_Y, p_base, 96) + beta_0 unique(X1.tips.exp.mat[1, ]) -# plot.data.process.actual(Y.state = X1.tips.exp.mat[1, ], +# plot_data.process.actual(Y.state = X1.tips.exp.mat[1, ], # phylo = trees[["96"]], # params = list(shifts = list(edges = shifts_grid[["96_3"]]$edges, values = shifts_grid[["96_3"]]$values[1, ])), # adj.root = 0, @@ -418,7 +418,7 @@ W <- PhylogeneticEM:::compute_actualization_matrix_ultrametric(trees[["160"]], a vec_Y <- kronecker(T_tree[["160"]], diag(1, p_base, p_base)) %*% W %*% as.vector(Delta) X1.tips.exp.mat <- matrix(vec_Y, p_base, 160) + beta_0 unique(X1.tips.exp.mat[1, ]) -# plot.data.process.actual(Y.state = X1.tips.exp.mat[1, ], +# plot_data.process.actual(Y.state = X1.tips.exp.mat[1, ], # phylo = trees[["160"]], # params = list(shifts = list(edges = shifts_grid[["160_3"]]$edges, values = shifts_grid[["160_3"]]$values[1, ])), # adj.root = 0, @@ -448,7 +448,7 @@ W <- PhylogeneticEM:::compute_actualization_matrix_ultrametric(trees[["160"]], a vec_Y <- kronecker(T_tree[["160"]], diag(1, p_base, p_base)) %*% W %*% as.vector(Delta) X1.tips.exp.mat <- matrix(vec_Y, p_base, 160) + beta_0 unique(X1.tips.exp.mat[1, ]) -# plot.data.process.actual(Y.state = X1.tips.exp.mat[1, ], +# plot_data.process.actual(Y.state = X1.tips.exp.mat[1, ], # phylo = trees[["160"]], # params = list(shifts = list(edges = shifts_grid[["160_7"]]$edges, values = shifts_grid[["160_7"]]$values[1, ])), # adj.root = 0, @@ -483,7 +483,7 @@ W <- PhylogeneticEM:::compute_actualization_matrix_ultrametric(trees[["160"]], a vec_Y <- kronecker(T_tree[["160"]], diag(1, p_base, p_base)) %*% W %*% as.vector(Delta) X1.tips.exp.mat <- matrix(vec_Y, p_base, 160) + beta_0 unique(X1.tips.exp.mat[1, ]) -# plot.data.process.actual(Y.state = X1.tips.exp.mat[1, ], +# plot_data.process.actual(Y.state = X1.tips.exp.mat[1, ], # phylo = trees[["160"]], # params = list(shifts = list(edges = shifts_grid[["160_11"]]$edges, values = shifts_grid[["160_11"]]$values[1, ])), # adj.root = 0, @@ -523,7 +523,7 @@ W <- PhylogeneticEM:::compute_actualization_matrix_ultrametric(trees[["160"]], a vec_Y <- kronecker(T_tree[["160"]], diag(1, p_base, p_base)) %*% W %*% as.vector(Delta) X1.tips.exp.mat <- matrix(vec_Y, p_base, 160) + beta_0 unique(X1.tips.exp.mat[1, ]) -# plot.data.process.actual(Y.state = X1.tips.exp.mat[1, ], +# plot_data.process.actual(Y.state = X1.tips.exp.mat[1, ], # phylo = trees[["160"]], # params = list(shifts = list(edges = shifts_grid[["160_15"]]$edges, values = shifts_grid[["160_15"]]$values[1, ])), # adj.root = 0, @@ -549,7 +549,7 @@ W <- PhylogeneticEM:::compute_actualization_matrix_ultrametric(trees[["192"]], a vec_Y <- kronecker(T_tree[["192"]], diag(1, p_base, p_base)) %*% W %*% as.vector(Delta) X1.tips.exp.mat <- matrix(vec_Y, p_base, 192) + beta_0 unique(X1.tips.exp.mat[1, ]) -# plot.data.process.actual(Y.state = X1.tips.exp.mat[1, ], +# plot_data.process.actual(Y.state = X1.tips.exp.mat[1, ], # phylo = trees[["192"]], # params = list(shifts = list(edges = shifts_grid[["192_3"]]$edges, values = shifts_grid[["192_3"]]$values[1, ])), # adj.root = 0, @@ -575,7 +575,7 @@ W <- PhylogeneticEM:::compute_actualization_matrix_ultrametric(trees[["256"]], a vec_Y <- kronecker(T_tree[["256"]], diag(1, p_base, p_base)) %*% W %*% as.vector(Delta) X1.tips.exp.mat <- matrix(vec_Y, p_base, 256) + beta_0 unique(X1.tips.exp.mat[1, ]) -# plot.data.process.actual(Y.state = X1.tips.exp.mat[1, ], +# plot_data.process.actual(Y.state = X1.tips.exp.mat[1, ], # phylo = trees[["256"]], # params = list(shifts = list(edges = shifts_grid[["256_3"]]$edges, values = shifts_grid[["256_3"]]$values[1, ])), # adj.root = 0, diff --git a/simulations_study/multivariate_exploitation.R b/simulations_study/multivariate_exploitation.R index 9aeb1e4..e5efbd9 100644 --- a/simulations_study/multivariate_exploitation.R +++ b/simulations_study/multivariate_exploitation.R @@ -1275,7 +1275,7 @@ simest$sim$gamma simest$sim$shifts simest$sim$rd -plot.data.process.actual(Y.state = simest$sim$Y_data[1,], +plot_data.process.actual(Y.state = simest$sim$Y_data[1,], phylo = trees[[paste0(simest$sim$ntaxa)]], params = params_test_true, adj.root = 1.3, @@ -1286,7 +1286,7 @@ plot.data.process.actual(Y.state = simest$sim$Y_data[1,], bg_beta_0 = "azure2", no.margin = TRUE) -plot.data.process.actual(Y.state = simest$sim$Y_data[1,], +plot_data.process.actual(Y.state = simest$sim$Y_data[1,], phylo = trees[[paste0(simest$sim$ntaxa)]], params = params_test, adj.root = 1.3, @@ -1297,7 +1297,7 @@ plot.data.process.actual(Y.state = simest$sim$Y_data[1,], bg_beta_0 = "azure2", no.margin = TRUE) -plot.data.process.actual(Y.state = simest$sim$Y_data[1,], +plot_data.process.actual(Y.state = simest$sim$Y_data[1,], phylo = trees[[paste0(simest$sim$ntaxa)]], params = simest$res$alpha_max$params_estim$`3`, adj.root = 1.3, @@ -1336,7 +1336,7 @@ sapply(results_estim_EM$params_history, function(z) z$shifts$edges) simest$sim$log_likelihood.true simest$res$alpha_max$results_summary$log_likelihood -plot.data.process.actual(Y.state = simest$sim$Y_data[1,], +plot_data.process.actual(Y.state = simest$sim$Y_data[1,], phylo = trees[[paste0(simest$sim$ntaxa)]], params = results_estim_EM_3$params, adj.root = 1.3, diff --git a/simulations_study/penalty_calibration.R b/simulations_study/penalty_calibration.R index b422128..6909e44 100644 --- a/simulations_study/penalty_calibration.R +++ b/simulations_study/penalty_calibration.R @@ -59,7 +59,7 @@ compute_distance <- function(simest){ Sigma <- compute_variance_covariance.OU(times_shared = times_shared, distances_phylo = distances_phylo, params_old = params_simu) - Sigma_YY <- extract.variance_covariance(Sigma, what="YY") + Sigma_YY <- extract_variance_covariance(Sigma, what="YY") Sigma_YY_inv <- solve(Sigma_YY) ## Mean at the tips (true parameters) XX <- simulate_internal(phylo = tree, diff --git a/simulations_study/several_K_data_generation.R b/simulations_study/several_K_data_generation.R index 733ea4e..9c2a1f7 100644 --- a/simulations_study/several_K_data_generation.R +++ b/simulations_study/several_K_data_generation.R @@ -130,7 +130,7 @@ datasetsim <- function(alpha, gamma, K, ntaxa, n, grp) { Sigma <- compute_variance_covariance.OU(times_shared = times_shared[[paste0(ntaxa)]], distances_phylo = distances_phylo[[paste0(ntaxa)]], params_old = params) - Sigma_YY <- extract.variance_covariance(Sigma, what="YY") + Sigma_YY <- extract_variance_covariance(Sigma, what="YY") Sigma_YY_inv <- solve(Sigma_YY) # Difficulty mu_0 <- (sum(Sigma_YY_inv))^(-1) * sum(Sigma_YY_inv%*%sim$m_Y_data) diff --git a/simulations_study/several_K_data_generation_uncertainties.R b/simulations_study/several_K_data_generation_uncertainties.R index fba2c7e..b5a14d5 100644 --- a/simulations_study/several_K_data_generation_uncertainties.R +++ b/simulations_study/several_K_data_generation_uncertainties.R @@ -148,7 +148,7 @@ datasetsim <- function(alpha, gamma, K, s_noise, ntaxa, n, grp) { Sigma <- compute_variance_covariance.scOU(times_shared = times_shared[[paste0(ntaxa)]], distances_phylo = distances_phylo[[paste0(ntaxa)]], params_old = params) - Sigma_YY <- extract.variance_covariance(Sigma, what="YY") + Sigma_YY <- extract_variance_covariance(Sigma, what="YY") Sigma_YY_chol <- chol(Sigma_YY) Sigma_YY_chol_inv <- backsolve(Sigma_YY_chol, diag(ncol(Sigma_YY_chol))) Sigma_YY_inv <- Sigma_YY_chol_inv %*% t(Sigma_YY_chol_inv) diff --git a/simulations_study/several_K_exploitation.R b/simulations_study/several_K_exploitation.R index 9940489..52316e9 100644 --- a/simulations_study/several_K_exploitation.R +++ b/simulations_study/several_K_exploitation.R @@ -558,14 +558,14 @@ conf <- 17 pars <- simestimations_alpha_known_K_select[[conf]] ## Clustering of tips and edges -plot.data.process.actual(Y.state = pars$m_Y_estim, +plot_data.process.actual(Y.state = pars$m_Y_estim, phylo = trees[[paste(pars$ntaxa)]], params = list(shifts = pars$params_estim$shifts, optimal.value = pars$params_estim$optimal.value), adj = 2, automatic_colors = TRUE) pars <- simlist[[conf]] -plot.data.process.actual(Y.state = pars$m_Y_data, +plot_data.process.actual(Y.state = pars$m_Y_data, phylo = trees[[paste(pars$ntaxa)]], params = list(shifts = pars$shifts, optimal.value = beta_0), diff --git a/simulations_study/several_K_exploitation_FUN.R b/simulations_study/several_K_exploitation_FUN.R index 4ce6ffd..e8d2e47 100644 --- a/simulations_study/several_K_exploitation_FUN.R +++ b/simulations_study/several_K_exploitation_FUN.R @@ -669,14 +669,14 @@ conf <- 17 pars <- simestimations_alpha_known_K_select[[conf]] ## Clustering of tips and edges -plot.data.process.actual(Y.state = pars$m_Y_estim, +plot_data.process.actual(Y.state = pars$m_Y_estim, phylo = trees[[paste(pars$ntaxa)]], params = list(shifts = pars$params_estim$shifts, optimal.value = pars$params_estim$optimal.value), adj = 2, automatic_colors = TRUE) pars <- simlist[[conf]] -plot.data.process.actual(Y.state = pars$m_Y_data, +plot_data.process.actual(Y.state = pars$m_Y_data, phylo = trees[[paste(pars$ntaxa)]], params = list(shifts = pars$shifts, optimal.value = beta_0), diff --git a/simulations_study/simulation_exploitation_bayou_design.R b/simulations_study/simulation_exploitation_bayou_design.R index de79555..d51eaa2 100644 --- a/simulations_study/simulation_exploitation_bayou_design.R +++ b/simulations_study/simulation_exploitation_bayou_design.R @@ -42,7 +42,7 @@ compute_difficulty <- function(simest){ Sigma <- compute_variance_covariance.OU(times_shared = times_shared, distances_phylo = distances_phylo, params_old = params_simu) - Sigma_YY <- extract.variance_covariance(Sigma, what="YY") + Sigma_YY <- extract_variance_covariance(Sigma, what="YY") Sigma_YY_inv <- solve(Sigma_YY) ## Mean at the tips XX <- simulate_internal(phylo = tree, diff --git a/simulations_study/simulation_exploitation_bayou_design_alpha_known.R b/simulations_study/simulation_exploitation_bayou_design_alpha_known.R index faf8ef3..cafb000 100644 --- a/simulations_study/simulation_exploitation_bayou_design_alpha_known.R +++ b/simulations_study/simulation_exploitation_bayou_design_alpha_known.R @@ -39,7 +39,7 @@ compute_difficulty <- function(simest){ Sigma <- compute_variance_covariance.OU(times_shared = times_shared, distances_phylo = distances_phylo, params_old = params_simu) - Sigma_YY <- extract.variance_covariance(Sigma, what="YY") + Sigma_YY <- extract_variance_covariance(Sigma, what="YY") Sigma_YY_inv <- solve(Sigma_YY) ## Mean at the tips XX <- simulate_internal(phylo = tree, diff --git a/tests/testthat/test-core-log_likelihood.R b/tests/testthat/test-core-log_likelihood.R index b392fa0..4f9f773 100644 --- a/tests/testthat/test-core-log_likelihood.R +++ b/tests/testthat/test-core-log_likelihood.R @@ -75,7 +75,7 @@ test_that("log-likelihood - scOU - random root", { params <- params_process.PhyloEM(res_new, K = 1, alpha = "max") params$selection.strength <- unique(diag(params$selection.strength)) C <- compute_tree_correlations_matrix.scOU(times_shared, distances_phylo, params) - C <- 1/(2*selection.strength) * extract.variance_covariance(C, what="YY", + C <- 1/(2*selection.strength) * extract_variance_covariance(C, what="YY", masque_data = c(rep(TRUE, ntaxa), rep(FALSE, dim(C)[1] - ntaxa))) C_inv <- solve(C) @@ -175,7 +175,7 @@ test_that("log-likelihood - BM - p=1 - fixed root", { params <- params_process.PhyloEM(res_new, K = 1, alpha = "max") params$selection.strength <- unique(diag(params$selection.strength)) C <- compute_tree_correlations_matrix.BM(times_shared, params) - C <- extract.variance_covariance(C, what="YY", + C <- extract_variance_covariance(C, what="YY", masque_data = c(rep(TRUE, ntaxa), rep(FALSE, dim(C)[1] - ntaxa))) C_inv <- solve(C) From 06a9b309736ef61801f0e2505371aa52ff3832fb Mon Sep 17 00:00:00 2001 From: Paul Bastide Date: Wed, 31 Jan 2024 16:37:30 +0100 Subject: [PATCH 08/10] Citation style --- inst/CITATION | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/inst/CITATION b/inst/CITATION index 54f2a63..2b22232 100644 --- a/inst/CITATION +++ b/inst/CITATION @@ -1,8 +1,8 @@ citHeader("To cite PhylogeneticEM in publications use:") -citEntry(entry = "Article", +bibentry(bibtype = "Article", title = "{Detection of adaptive shifts on phylogenies by using shifted stochastic processes on a tree}", - author = personList(as.person("Paul Bastide"), + author = c(as.person("Paul Bastide"), as.person("Mahendra Mariadassou"), as.person("Stephane Robin")), journal = "Journal of the Royal Statistical Society: Series B (Statistical Methodology)", From b39ab3cf01336f9f109e154869b46c30b400a853 Mon Sep 17 00:00:00 2001 From: Paul Bastide Date: Wed, 31 Jan 2024 16:39:58 +0100 Subject: [PATCH 09/10] No parallel computing in vignettes. --- vignettes/monkeys.Rmd | 119 ++++++++++++----------------------------- vignettes/tutorial.Rmd | 32 ++++------- 2 files changed, 43 insertions(+), 108 deletions(-) diff --git a/vignettes/monkeys.Rmd b/vignettes/monkeys.Rmd index 8cf6bcd..56b8f55 100644 --- a/vignettes/monkeys.Rmd +++ b/vignettes/monkeys.Rmd @@ -7,7 +7,6 @@ vignette: > %\VignetteIndexEntry{New World Monkeys} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} - %\VignetteDepends{doParallel} --- ## Dataset @@ -49,26 +48,14 @@ res ``` ```{r label="Fit_EM_int", echo=FALSE, warning=FALSE, message=FALSE} -if (!requireNamespace("doParallel", quietly = TRUE)) { - res <- PhyloEM(phylo = monkeys$phy, - Y_data = monkeys$dat, - process = "scOU", ## scalar OU model - random.root = TRUE, ## Root is stationary - stationary.root = TRUE, - nbr_alpha = 4, ## Number of alpha values tested (should be raised) - K_max = 10, ## Maximal number of shifts - parallel_alpha = FALSE) ## This can be set to TRUE for -} else { - res <- PhyloEM(phylo = monkeys$phy, - Y_data = monkeys$dat, - process = "scOU", ## scalar OU model - random.root = TRUE, ## Root is stationary - stationary.root = TRUE, - nbr_alpha = 4, ## Number of alpha values tested (should be raised) - K_max = 10, ## Maximal number of shifts - parallel_alpha = TRUE, ## This can be set to TRUE for - Ncores = 2) ## parallel computations -} +res <- PhyloEM(phylo = monkeys$phy, + Y_data = monkeys$dat, + process = "scOU", ## scalar OU model + random.root = TRUE, ## Root is stationary + stationary.root = TRUE, + nbr_alpha = 4, ## Number of alpha values tested (should be raised) + K_max = 10, ## Maximal number of shifts + parallel_alpha = FALSE) ## This can be set to TRUE for res ``` @@ -153,28 +140,15 @@ rot <- matrix(c(cos(theta), -sin(theta), sin(theta), cos(theta)), Yrot <- t(rot) %*% monkeys$dat rownames(Yrot) <- rownames(monkeys$dat) -if (!requireNamespace("doParallel", quietly = TRUE)) { - # PhyloEM on rotated data - res_rot <- PhyloEM(phylo = monkeys$phy, - Y_data = Yrot, ## rotated data - process = "scOU", - random.root = TRUE, - stationary.root = TRUE, - nbr_alpha = 4, - K_max = 10, - parallel_alpha = FALSE) -} else { - # PhyloEM on rotated data - res_rot <- PhyloEM(phylo = monkeys$phy, - Y_data = Yrot, ## rotated data - process = "scOU", - random.root = TRUE, - stationary.root = TRUE, - nbr_alpha = 4, - K_max = 10, - parallel_alpha = TRUE, - Ncores = 2) -} +# PhyloEM on rotated data +res_rot <- PhyloEM(phylo = monkeys$phy, + Y_data = Yrot, ## rotated data + process = "scOU", + random.root = TRUE, + stationary.root = TRUE, + nbr_alpha = 4, + K_max = 10, + parallel_alpha = FALSE) ``` ```{r, fig.show='hold', fig.height=4, fig.width=7, warning=FALSE} @@ -262,28 +236,15 @@ rot2 <- matrix(c(cos(theta), -sin(theta), sin(theta), cos(theta)), Yrot2 <- t(rot2) %*% monkeys$dat rownames(Yrot2) <- rownames(monkeys$dat) -if (!requireNamespace("doParallel", quietly = TRUE)) { - # PhyloEM on rotated data - res_rot2 <- PhyloEM(phylo = monkeys$phy, - Y_data = Yrot2, - process = "scOU", - random.root = TRUE, - stationary.root = TRUE, - nbr_alpha = 2, ## Note that this can also be different - K_max = 10, - parallel_alpha = FALSE) -} else { - # PhyloEM on rotated data - res_rot2 <- PhyloEM(phylo = monkeys$phy, - Y_data = Yrot2, - process = "scOU", - random.root = TRUE, - stationary.root = TRUE, - nbr_alpha = 2, ## Note that this can also be different - K_max = 10, - parallel_alpha = TRUE, - Ncores = 2) -} +# PhyloEM on rotated data +res_rot2 <- PhyloEM(phylo = monkeys$phy, + Y_data = Yrot2, + process = "scOU", + random.root = TRUE, + stationary.root = TRUE, + nbr_alpha = 2, ## Note that this can also be different + K_max = 10, + parallel_alpha = FALSE) # Merge res_merge2 <- merge_rotations(res, res_rot, res_rot2) ``` @@ -325,26 +286,14 @@ res_neg <- PhyloEM(phylo = monkeys$phy, ``` ```{r label="Fit_EM_negative_int", echo=FALSE, warning=FALSE, message=FALSE} -if (!requireNamespace("doParallel", quietly = TRUE)) { - res_neg <- PhyloEM(phylo = monkeys$phy, - Y_data = monkeys$dat, - process = "scOU", - random.root = FALSE, ## root needs to be fixed - K_max = 10, - parallel_alpha = FALSE, - nbr_alpha = 4, ## 2 negative, 2 positive (should be more) - allow_negative = TRUE) ## allow negative alpha in the grid -} else { - res_neg <- PhyloEM(phylo = monkeys$phy, - Y_data = monkeys$dat, - process = "scOU", - random.root = FALSE, ## root needs to be fixed - K_max = 10, - parallel_alpha = TRUE, - Ncores = 2, - nbr_alpha = 4, ## 2 negative, 2 positive (should be more) - allow_negative = TRUE) ## allow negative alpha in the grid -} +res_neg <- PhyloEM(phylo = monkeys$phy, + Y_data = monkeys$dat, + process = "scOU", + random.root = FALSE, ## root needs to be fixed + K_max = 10, + parallel_alpha = FALSE, + nbr_alpha = 4, ## 2 negative, 2 positive (should be more) + allow_negative = TRUE) ## allow negative alpha in the grid ``` We can extract the parameter of the fit with no shift, that is equivalent to an EB model: diff --git a/vignettes/tutorial.Rmd b/vignettes/tutorial.Rmd index b79bb4a..5ab9e2f 100644 --- a/vignettes/tutorial.Rmd +++ b/vignettes/tutorial.Rmd @@ -8,7 +8,6 @@ vignette: > %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} %\VignetteDepends{TreeSim} - %\VignetteDepends{doParallel} --- ## Introduction @@ -131,28 +130,15 @@ res ```{r label="Fit_EM_int", echo=FALSE, warning = FALSE} ## Grid on alpha alpha_grid <- c(1, 3) -if (!requireNamespace("doParallel", quietly = TRUE)) { - ## Run algorithm - res <- PhyloEM(phylo = tree, - Y_data = data, - process = "scOU", ## scalar OU model - random.root = TRUE, ## Root is stationary (true model) - stationary.root = TRUE, - alpha = alpha_grid, ## On a grid of alpha - K_max = 10, ## Maximal number of shifts - parallel_alpha = FALSE) ## This can be set to TRUE for -} else { - ## Run algorithm - res <- PhyloEM(phylo = tree, - Y_data = data, - process = "scOU", ## scalar OU model - random.root = TRUE, ## Root is stationary (true model) - stationary.root = TRUE, - alpha = alpha_grid, ## On a grid of alpha - K_max = 10, ## Maximal number of shifts - parallel_alpha = TRUE, ## This can be set to TRUE for - Ncores = 2) ## parallel computations -} +## Run algorithm +res <- PhyloEM(phylo = tree, + Y_data = data, + process = "scOU", ## scalar OU model + random.root = TRUE, ## Root is stationary (true model) + stationary.root = TRUE, + alpha = alpha_grid, ## On a grid of alpha + K_max = 10, ## Maximal number of shifts + parallel_alpha = FALSE) ## This can be set to TRUE for res ``` From 09f8f12a38a245c1dbc8bfbfd69f3050593c6682 Mon Sep 17 00:00:00 2001 From: Paul Bastide Date: Wed, 31 Jan 2024 17:02:27 +0100 Subject: [PATCH 10/10] Getting ready for CRAN --- DESCRIPTION | 2 +- NEWS.md | 2 +- cran-comments.md | 26 +++----------------------- 3 files changed, 5 insertions(+), 25 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index bb31afe..2281950 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: PhylogeneticEM Title: Automatic Shift Detection using a Phylogenetic EM -Version: 1.6.0.9000 +Version: 1.7.0 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 714e5d8..0cf00b9 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,4 +1,4 @@ -# PhylogeneticEM 1.6.0.9000 +# PhylogeneticEM 1.7.0 * Bug fixes: * Fix documentation notes on CRAN checks. * New functions and features: diff --git a/cran-comments.md b/cran-comments.md index 3fd4cb9..1bba000 100644 --- a/cran-comments.md +++ b/cran-comments.md @@ -1,5 +1,5 @@ ## Test environments -* local: macOS 10.15, R 4.2.1 +* local: macOS 14.2.1, R 4.3.2 * GitHub Actions: * macOS-latest: release * windows-latest: release @@ -7,27 +7,7 @@ * win-builder: devel and release ## R CMD check results -There were no ERRORs or WARNINGs. -There were two NOTEs: -``` -checking installed package size ... NOTE - installed size is 5.0Mb - sub-directories of 1Mb or more: - libs 3.5Mb -``` -This package uses `Rcpp` and `RcppArmadillo`. - -``` -* checking CRAN incoming feasibility ... [18s] NOTE -Maintainer: 'Paul Bastide ' - -Found the following (possibly) invalid URLs: - URL: https://doi.org/10.1111/rssb.12206 - From: inst/CITATION - Status: 503 - Message: Service Unavailable -``` -This doi and this address should be valid. +There were no ERROR, WARNING or NOTE. ## Downstream dependencies -There are currently no downstream dependencies for this package. \ No newline at end of file +There are currently no downstream dependencies for this package.