diff --git a/DESCRIPTION b/DESCRIPTION index 83183c6..c055a0e 100755 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -2,7 +2,7 @@ Package: secsse Type: Package Title: Several Examined and Concealed States-Dependent Speciation and Extinction -Version: 3.6.0 +Version: 3.6.1 License: GPL (>= 3) | file LICENSE Authors@R: c( person(given = "Leonel", @@ -50,7 +50,8 @@ Imports: ggplot2, tibble, rlang, - treestats + treestats, + pracma Suggests: diversitree, phytools, diff --git a/NEWS.md b/NEWS.md index 94a3ee4..2901193 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,10 +1,18 @@ +# 3.6.1 +- added support for having a different number of concealed states compared to +the number of observed states when creating a q_matrix using `q_doubletrans`. + # 3.6.0 - changed the default integration method to "odeint::runge_kutta_cash_karp54", although a bit slower, this method is more numerically stable. - fixed lack of sorting of numeric traits in the function `q_doubletrans` -- fixed an error in preparing the state matrices when NAs were present in the traits -- ML and LL functions now optionally also return the root state, which in turn can be used to in `secsse_sim` as a starting point at the root / crown of the tree. -- updated simulations to sample species using a binary search, instead of using stochastic acceptance +- fixed an error in preparing the state matrices when NAs were present in the +traits +- ML and LL functions now optionally also return the root state, which in turn +can be used to in `secsse_sim` as a starting point at the root / crown of the +tree. +- updated simulations to sample species using a binary search, instead of using +stochastic acceptance - better documentation of output of ML and LL functions. # 3.5.0 Version 3.5.0 uses a separate calculation for 1 - E, e.g. one minus the local diff --git a/R/default_params_doc.R b/R/default_params_doc.R index 6c4f616..ebb9784 100644 --- a/R/default_params_doc.R +++ b/R/default_params_doc.R @@ -33,7 +33,8 @@ #' @param cond condition on the existence of a node root: `"maddison_cond"`, #' `"proper_cond"` (default). For details, see vignette. #' @param root_state_weight the method to weigh the states: -#' `"maddison_weights"`, `"proper_weights"` (default) or `"equal_weights"`. +#' `"maddison_weights"`, `"proper_weights"` (default), `"equal_weights"'. +#' or `"stationary_weights"` #' It can also be specified for the root state: the vector `c(1, 0, 0)` #' indicates state 1 was the root state. When #' using a `multiPhylo` object, root_state_weight should be list where each diff --git a/R/seccse_multi_loglik.R b/R/seccse_multi_loglik.R index 114b895..4901426 100644 --- a/R/seccse_multi_loglik.R +++ b/R/seccse_multi_loglik.R @@ -45,31 +45,31 @@ multi_loglik <- function(parameter, } if (length(phy[[i]]$tip.label) == 1) { - local_res <- secsse::secsse_single_branch_loglik(parameter = parameter, - phy = phy[[i]], - traits = traits[[i]], - num_concealed_states = - num_concealed_states, - cond = cond, - root_state_weight = - focal_root_state_weight, - sampling_fraction = - focal_sampling_fraction, - setting_calculation = - focal_setting_calculation, - see_ancestral_states = FALSE, - loglik_penalty = loglik_penalty, - is_complete_tree = - is_complete_tree, - take_into_account_root_edge = - take_into_account_root_edge, - num_threads = num_threads, - atol = atol, - rtol = rtol, - method = method, - display_warning = display_warning, - use_normalization = use_normalization, - return_root_state = return_root_state) + local_res <- secsse_single_branch_loglik(parameter = parameter, + phy = phy[[i]], + traits = traits[[i]], + num_concealed_states = + num_concealed_states, + cond = cond, + root_state_weight = + focal_root_state_weight, + sampling_fraction = + focal_sampling_fraction, + setting_calculation = + focal_setting_calculation, + see_ancestral_states = FALSE, + loglik_penalty = loglik_penalty, + is_complete_tree = + is_complete_tree, + take_into_account_root_edge = + take_into_account_root_edge, + num_threads = num_threads, + atol = atol, + rtol = rtol, + method = method, + display_warning = display_warning, + use_normalization = use_normalization, + return_root_state = return_root_state) if (return_root_state) { root_states[[i]] <- local_res$root_state } @@ -107,6 +107,8 @@ multi_loglik <- function(parameter, ll <- do.call(sum, res) + + if (return_root_state) { return(list(LL = ll, root_state = root_states)) diff --git a/R/secsse_loglik.R b/R/secsse_loglik.R index 85cb0c9..ba99636 100644 --- a/R/secsse_loglik.R +++ b/R/secsse_loglik.R @@ -54,32 +54,32 @@ master_loglik <- function(parameter, } if (length(phy$tip.label) == 1) { - return(secsse::secsse_single_branch_loglik( - parameter = parameter, - phy = phy, - traits = traits, - num_concealed_states = - num_concealed_states, - cond = cond, - root_state_weight = - root_state_weight, - sampling_fraction = - sampling_fraction, - setting_calculation = - setting_calculation, - see_ancestral_states = see_ancestral_states, - loglik_penalty = loglik_penalty, - is_complete_tree = - is_complete_tree, - take_into_account_root_edge = - take_into_account_root_edge, - num_threads = num_threads, - atol = atol, - rtol = rtol, - method = method, - display_warning = display_warning, - use_normalization = use_normalization, - return_root_state = return_root_state)) + return(secsse_single_branch_loglik( + parameter = parameter, + phy = phy, + traits = traits, + num_concealed_states = + num_concealed_states, + cond = cond, + root_state_weight = + root_state_weight, + sampling_fraction = + sampling_fraction, + setting_calculation = + setting_calculation, + see_ancestral_states = see_ancestral_states, + loglik_penalty = loglik_penalty, + is_complete_tree = + is_complete_tree, + take_into_account_root_edge = + take_into_account_root_edge, + num_threads = num_threads, + atol = atol, + rtol = rtol, + method = method, + display_warning = display_warning, + use_normalization = use_normalization, + return_root_state = return_root_state)) } lambdas <- parameter[[1]] @@ -113,7 +113,6 @@ master_loglik <- function(parameter, states <- setting_calculation$states forTime <- setting_calculation$forTime ances <- setting_calculation$ances - d <- ncol(states) / 3 # with a complete tree, we need to re-calculate the states every time we @@ -193,7 +192,8 @@ master_loglik <- function(parameter, lambdas, nodeM, d, - is_cla = using_cla) + is_cla = using_cla, + Q = q_matrix) if (is_complete_tree) { nodeM <- update_complete_tree(phy, @@ -325,25 +325,26 @@ secsse_loglik <- function(parameter, display_warning = TRUE, use_normalization = TRUE, return_root_state = FALSE) { - master_loglik(parameter = parameter, - phy = phy, - traits = traits, - num_concealed_states = num_concealed_states, - cond = cond, - root_state_weight = root_state_weight, - sampling_fraction = sampling_fraction, - setting_calculation = setting_calculation, - see_ancestral_states = see_ancestral_states, - loglik_penalty = loglik_penalty, - is_complete_tree = is_complete_tree, - take_into_account_root_edge = take_into_account_root_edge, - num_threads = num_threads, - atol = atol, - rtol = rtol, - method = method, - display_warning = display_warning, - use_normalization = use_normalization, - return_root_state = return_root_state) + ll <- master_loglik(parameter = parameter, + phy = phy, + traits = traits, + num_concealed_states = num_concealed_states, + cond = cond, + root_state_weight = root_state_weight, + sampling_fraction = sampling_fraction, + setting_calculation = setting_calculation, + see_ancestral_states = see_ancestral_states, + loglik_penalty = loglik_penalty, + is_complete_tree = is_complete_tree, + take_into_account_root_edge = take_into_account_root_edge, + num_threads = num_threads, + atol = atol, + rtol = rtol, + method = method, + display_warning = display_warning, + use_normalization = use_normalization, + return_root_state = return_root_state) + return(ll) } diff --git a/R/secsse_ml.R b/R/secsse_ml.R index b6cd4e8..f8cb25a 100644 --- a/R/secsse_ml.R +++ b/R/secsse_ml.R @@ -29,6 +29,7 @@ master_ml <- function(phy, use_normalization = TRUE, return_root_state = FALSE, see_ancestral_states = FALSE) { + structure_func <- NULL if (!is.null(functions_defining_params)) { structure_func <- set_and_check_structure_func(idparsfuncdefpar, @@ -81,16 +82,6 @@ master_ml <- function(phy, num_modeled_traits <- length(idparslist[[1]]) / num_concealed_states if (!is.list(traits)) { - - if (length(phy$tip.label) == 1) { - fake_phy <- ape::rphylo(n = 2, birth = 1, death = 0) - fake_phy$edge.length[1:2] <- phy$edge.length[1] - old_phy <- phy - old_traits <- traits - phy <- fake_phy - traits <- c(traits, traits) - } - setting_calculation <- build_initStates_time(phy, traits, num_concealed_states, @@ -100,17 +91,11 @@ master_ml <- function(phy, num_modeled_traits, traitStates = get_trait_states(idparslist, - num_concealed_states, FALSE)) - if (length(phy$tip.label) == 1) { - setting_calculation$states <- setting_calculation$states[-2, ] - setting_calculation$forTime <- setting_calculation$forTime[-2, ] - phy <- old_phy - traits <- old_traits - } + num_concealed_states, verbose)) + } else { setting_calculation <- list() for (i in 1:length(phy)) { - input_phy <- phy[[i]] input_traits <- traits[[i]] @@ -119,14 +104,7 @@ master_ml <- function(phy, } else { input_sampling_fraction <- sampling_fraction } - - if (length(input_phy$tip.label) == 1) { - fake_phy <- ape::rphylo(n = 2, birth = 1, death = 0) - fake_phy$edge.length[1:2] <- input_phy$edge.length[1] - input_phy <- fake_phy - input_traits <- c(input_traits, input_traits) - } - + setting_calculation[[i]] <- build_initStates_time(phy = input_phy, traits = input_traits, num_concealed_states = @@ -141,7 +119,7 @@ master_ml <- function(phy, first_time = FALSE, traitStates = get_trait_states(idparslist, - num_concealed_states, FALSE)) + num_concealed_states, verbose)) } } diff --git a/R/secsse_prep.R b/R/secsse_prep.R index fcbaeb5..1564adf 100644 --- a/R/secsse_prep.R +++ b/R/secsse_prep.R @@ -145,9 +145,7 @@ create_q_matrix <- function(state_names, diff.conceal = FALSE) { total_num_states <- length(state_names) - if (total_num_states != num_concealed_states) { - stop("number of concealed states has to be equal to the number of observed states") - } + trans_matrix <- matrix(0, ncol = total_num_states, nrow = total_num_states) @@ -161,9 +159,10 @@ create_q_matrix <- function(state_names, diag(trans_matrix) <- NA - trans_matrix <- secsse::q_doubletrans(traits = state_names, - masterBlock = trans_matrix, - diff.conceal = diff.conceal) + trans_matrix <- q_doubletrans(traits = state_names, + masterBlock = trans_matrix, + diff.conceal = diff.conceal, + num_concealed_states = num_concealed_states) all_state_names <- get_state_names(state_names, num_concealed_states) colnames(trans_matrix) <- all_state_names @@ -190,9 +189,9 @@ expand_q_matrix <- function(q_matrix, traits <- get_state_names(names(q_matrix), num_concealed_states) - return(secsse::q_doubletrans(traits = traits, - masterBlock = q_matrix, - diff.conceal = diff.conceal)) + return(q_doubletrans(traits = traits, + masterBlock = q_matrix, + diff.conceal = diff.conceal)) } #' Helper function to create a default `shift_matrix` list @@ -310,36 +309,16 @@ create_mu_vector <- function(state_names, #' @keywords internal replace_matrix <- function(focal_matrix, - params, - is_lambda = FALSE) { - if (is_lambda) { - entries <- table(focal_matrix) - entries <- entries[names(entries) != 0] - for (i in seq_len(nrow(focal_matrix))) { - for (j in seq_len(ncol(focal_matrix))) { - if (focal_matrix[i, j] != 0 && !is.na(focal_matrix[i, j])) { - index <- focal_matrix[i, j] - - new_val <- params[index] - if (index %in% names(entries)) { - mult <- 1 / entries[names(entries) == index] - new_val <- new_val * mult - } - - focal_matrix[i, j] <- new_val - } - } - } - } else { - for (i in seq_len(nrow(focal_matrix))) { - for (j in seq_len(ncol(focal_matrix))) { - if (focal_matrix[i, j] != 0 && !is.na(focal_matrix[i, j])) { - index <- focal_matrix[i, j] - focal_matrix[i, j] <- params[index] - } + params) { + for (i in seq_len(nrow(focal_matrix))) { + for (j in seq_len(ncol(focal_matrix))) { + if (focal_matrix[i, j] != 0 && !is.na(focal_matrix[i, j])) { + index <- focal_matrix[i, j] + focal_matrix[i, j] <- params[index] } } } + return(focal_matrix) } @@ -353,7 +332,7 @@ fill_in <- function(object, params) { if (is.list(object)) { # lambda matrix for (k in seq_along(object)) { - object[[k]] <- replace_matrix(object[[k]], params, is_lambda = FALSE) + object[[k]] <- replace_matrix(object[[k]], params) } } else if (is.matrix(object)) { object <- replace_matrix(object, params) diff --git a/R/secsse_sim.R b/R/secsse_sim.R index 5517fda..9e66702 100644 --- a/R/secsse_sim.R +++ b/R/secsse_sim.R @@ -42,16 +42,16 @@ secsse_sim <- function(lambdas, if (is.matrix(lambdas)) { # need to be converted lambdas <- prepare_full_lambdas(names(mus), - num_concealed_states = num_concealed_states, - lambdas) + num_concealed_states = num_concealed_states, + lambdas) } - + if (length(lambdas) != length(mus)) { stop("Every state must have a single rate of speciation and extinction") } - + diag(qs) <- 0 - + if (is.null(init_state_probs)) { num_init_states <- length(mus) init_state_probs <- rep(1 / num_init_states, num_init_states) @@ -73,9 +73,9 @@ secsse_sim <- function(lambdas, init_state_probs <- init_state_probs / sum(init_state_probs) } - + if (is.null(seed)) seed <- -1 - + condition_vec <- vector() if (length(conditioning) > 1) { condition_vec <- conditioning @@ -91,26 +91,27 @@ secsse_sim <- function(lambdas, condition_vec <- -1 + indices } - while(TRUE) { + total_tries <- 0 + while(total_tries < max_tries) { res <- generate_phy(mus, - lambdas, - qs, - crown_age, - max_spec, - max_species_extant, - min_spec, - init_state_probs, - conditioning, - num_concealed_states, - non_extinction, - verbose, - max_tries, - seed, - condition_vec, - tree_size_hist, - start_at_crown, - drop_extinct) - + lambdas, + qs, + crown_age, + max_spec, + max_species_extant, + min_spec, + init_state_probs, + conditioning, + num_concealed_states, + non_extinction, + verbose, + max_tries, + seed, + condition_vec, + tree_size_hist, + start_at_crown, + drop_extinct) + total_tries <- total_tries + res$total_tries if (res$status == "not enough tries") { return(res) } @@ -118,7 +119,7 @@ secsse_sim <- function(lambdas, if (length(res$phy$tip.label) != length(res$obs_traits)) { message("something went wrong, please report this to the package maintainer") } - + if (res$status == "success" || res$status == "single_species_tree") { if (sum(sampling_fraction) == length(sampling_fraction) || is.null(sampling_fraction)) { @@ -141,18 +142,25 @@ secsse_sim <- function(lambdas, } # now we need to remove all the tips if (length(tips_to_remove) > 0) { - res$obs_traits <- res$obs_traits[-tips_to_remove] - res$true_traits <- res$true_traits[-tips_to_remove] - res$phy <- ape::drop.tip(res$phy, tips_to_remove) + if (length(tips_to_remove) == length(res$phy$tip.label)) { + res$phy <- NULL + res$obs_traits <- NULL + res$true_traits <- NULL + } else { + res$obs_traits <- res$obs_traits[-tips_to_remove] + res$true_traits <- res$true_traits[-tips_to_remove] + + res$phy <- ape::drop.tip(res$phy, tips_to_remove) + } } - + if (is.null(res$phy)) { warning("sampling removed all tips, trying again with new seed") seed = sample(1:1e9, 1) } else { return(res) } - } + } } return(res) } @@ -197,17 +205,16 @@ generate_phy <- function(mus, if (length(res) < 1) { # this happens upon a throw return(list(phy = "ds", traits = 0, - status = "error")) + status = "error", + total_tries = max_tries)) } - - Ltable <- res$ltable out_hist <- 0 if (tree_size_hist == TRUE) out_hist <- res$hist_tree_size - if (sum(res$tracker) >= max_tries) { + if (sum(res$tracker) > max_tries) { warning("Couldn't simulate a tree in enough tries, try increasing max_tries") @@ -218,7 +225,8 @@ generate_phy <- function(mus, conditioning = res$tracker[4], small = res$tracker[6], size_hist = out_hist, - status = "not enough tries")) + status = "not enough tries", + total_tries = sum(res$tracker))) } if (start_at_crown == FALSE && sum(Ltable[, 4] == -1) == 1) { @@ -248,7 +256,8 @@ generate_phy <- function(mus, conditioning = res$tracker[4], small = res$tracker[6], size_hist = out_hist, - status = "single_species_tree")) + status = "single_species_tree", + total_tries = sum(res$tracker))) } else if (sum(Ltable[, 4] == -1) < 2) { @@ -260,7 +269,8 @@ generate_phy <- function(mus, conditioning = res$tracker[4], small = res$tracker[6], size_hist = out_hist, - status = "extinction")) + status = "extinction", + total_tries = sum(res$tracker))) } @@ -310,7 +320,8 @@ generate_phy <- function(mus, conditioning = res$tracker[4], small = res$tracker[6], size_hist = out_hist, - status = "success")) + status = "success", + total_tries = sum(res$tracker))) } else { warning("simulation did not meet minimal requirements") return(list(phy = "ds", @@ -320,6 +331,7 @@ generate_phy <- function(mus, conditioning = res$tracker[4], small = res$tracker[6], size_hist = out_hist, - status = "requirements not met")) + status = "requirements not met", + total_tries = sum(res$tracker))) } } \ No newline at end of file diff --git a/R/secsse_single_branch.R b/R/secsse_single_branch.R index 30c943a..9dd57b3 100644 --- a/R/secsse_single_branch.R +++ b/R/secsse_single_branch.R @@ -39,13 +39,9 @@ secsse_single_branch_loglik <- function(parameter, if (is.null(setting_calculation)) { check_root_state_weight(root_state_weight, traits) - - # make fake phy - fake_phy <- ape::rphylo(n = 2, birth = 1, death = 0) - fake_phy$edge.length[1:2] <- phy$edge.length[1] - - setting_calculation <- build_initStates_time(fake_phy, - c(traits, traits), + + setting_calculation <- build_initStates_time(phy, + traits, num_concealed_states, sampling_fraction, is_complete_tree, @@ -55,14 +51,12 @@ secsse_single_branch_loglik <- function(parameter, } states <- setting_calculation$states - states <- states[-2, ] forTime <- setting_calculation$forTime - forTime <- forTime[-2, ] - - d <- ncol(states) / 3 + + d <- length(states) / 3 if (!is.null(phy$root.edge)) { - forTime[3] <- forTime[3] + phy$root.edge + forTime[2] <- forTime[2] + phy$root.edge } RcppParallel::setThreadOptions(numThreads = num_threads) @@ -71,8 +65,8 @@ secsse_single_branch_loglik <- function(parameter, if (return_root_state) return_states = TRUE calcul <- calc_ll_single_branch_cpp(rhs = if (using_cla) "ode_cla" else "ode_standard", - states = states[1, ], - forTime = c(0, forTime[3]), + states = states, + forTime = forTime, lambdas = lambdas, mus = mus, Q = q_matrix, @@ -97,7 +91,8 @@ secsse_single_branch_loglik <- function(parameter, lambdas, nodeM, d, - is_cla = using_cla) + is_cla = using_cla, + Q = q_matrix) mergeBranch2 <- condition(cond, mergeBranch, weight_states, diff --git a/R/secsse_utils.R b/R/secsse_utils.R index cf36b7e..a566858 100755 --- a/R/secsse_utils.R +++ b/R/secsse_utils.R @@ -58,37 +58,6 @@ id_paramPos <- function(traits, num_concealed_states) { #noLint return(idparslist) } -#' @keywords internal -create_q_matrix_int <- function(masterBlock, - concealnewQMatr, - ntraits, - diff.conceal) { - Q <- NULL - for (i in 1:ntraits) { - Qrow <- NULL - for (ii in 1:ntraits) { - entry <- masterBlock[i, ii] - if (is.na(entry)) { - Qrow <- cbind(Qrow, masterBlock) - } else { - if (diff.conceal == TRUE) { - entry <- concealnewQMatr[i, ii] - } - - outDiagBlock <- matrix(0, - ncol = ntraits, - nrow = ntraits, - byrow = TRUE) - diag(outDiagBlock) <- entry - Qrow <- cbind(Qrow, outDiagBlock) - } - } - Q <- rbind(Q, Qrow) - } - return(Q) -} - - #' @title Basic Qmatrix #' Sets a Q matrix where double transitions are not allowed #' @@ -96,90 +65,90 @@ create_q_matrix_int <- function(masterBlock, #' #' @return Q matrix that includes both examined and concealed states, it should #' be declared as the third element of idparslist. -#' @description This function expands the Q_matrix, but it does so assuming -#' that the number of concealed traits is equal to the number of examined -#' traits, if you have a different number, you should consider looking at -#' the function [expand_q_matrix()]. +#' @description This function expands the Q_matrix. If the number of concealed +#' states is not explicitly set by the user, it is assumed to be identical +#' to the number of observed states. #' @examples -#' traits <- sample(c(0,1,2), 45,replace = TRUE) #get some traits +#' traits <- sample(c(0, 1, 2), 45, replace = TRUE) # get some traits #' # For a three-state trait -#' masterBlock <- matrix(99,ncol = 3,nrow = 3,byrow = TRUE) +#' masterBlock <- matrix(99, ncol = 3, nrow = 3, byrow = TRUE) #' diag(masterBlock) <- NA -#' masterBlock[1,2] <- 6 -#' masterBlock[1,3] <- 7 -#' masterBlock[2,1] <- 8 -#' masterBlock[2,3] <- 9 -#' masterBlock[3,1] <- 10 -#' masterBlock[3,2] <- 11 -#' myQ <- q_doubletrans(traits,masterBlock,diff.conceal = FALSE) +#' masterBlock[1, 2] <- 6 +#' masterBlock[1, 3] <- 7 +#' masterBlock[2, 1] <- 8 +#' masterBlock[2, 3] <- 9 +#' masterBlock[3, 1] <- 10 +#' masterBlock[3, 2] <- 11 +#' myQ <- q_doubletrans(traits, masterBlock, diff.conceal = FALSE) #' # now, it can replace the Q matrix from id_paramPos #' num_concealed_states <- 3 #' param_posit <- id_paramPos(traits,num_concealed_states) #' param_posit[[3]] <- myQ #' @export -q_doubletrans <- function(traits, masterBlock, diff.conceal) { - if (diff.conceal == TRUE && - all(floor(masterBlock) == masterBlock, na.rm = TRUE) == FALSE) { - integersmasterBlock <- floor(masterBlock) - factorBlock <- signif(masterBlock - integersmasterBlock, digits = 2) - - factorstoExpand <- unique(sort(c(factorBlock))) - factorstoExpand <- factorstoExpand[factorstoExpand > 0] - newshareFac <- - (max(factorstoExpand * 10) + 1):(max(factorstoExpand * 10) + - length(factorstoExpand)) - newshareFac <- newshareFac / 10 - - for (iii in seq_along(newshareFac)) { - factorBlock[which(factorBlock == factorstoExpand[iii])] <- - newshareFac[iii] - } +q_doubletrans <- function(traits, + masterBlock, + diff.conceal, + num_concealed_states = NULL) { - ntraits <- length(sort(unique(traits))) - uniqParQ <- sort(unique(c(floor(masterBlock)))) - uniqParQ2 <- uniqParQ[which(uniqParQ > 0)] - concealnewQ <- (max(uniqParQ2) + 1):(max(uniqParQ2) + length(uniqParQ2)) - - for (iii in seq_along(concealnewQ)) { - integersmasterBlock[which(integersmasterBlock == uniqParQ2[iii])] <- - concealnewQ[iii] - } - concealnewQMatr <- integersmasterBlock + factorBlock + n_obs_traits <- length(sort(unique(traits))) + if (is.null(num_concealed_states)) num_concealed_states <- n_obs_traits + + Q <- matrix(data = 0, + nrow = n_obs_traits * num_concealed_states, + ncol = n_obs_traits * num_concealed_states) + + uniq_obs_rates <- sort(unique(c(masterBlock))) + uniq_obs_rates <- uniq_obs_rates[!is.na(uniq_obs_rates)] + uniq_obs_rates <- uniq_obs_rates[uniq_obs_rates > 0] + max_obs_rate <- max(uniq_obs_rates) + new_conceal_rates <- + (max_obs_rate + 1):(max_obs_rate + length(uniq_obs_rates)) + if (diff.conceal == FALSE) { + new_conceal_rates <- uniq_obs_rates + } - Q <- create_q_matrix_int(masterBlock, - concealnewQMatr, - ntraits, - diff.conceal) - } else { - ntraits <- length(sort(unique(traits))) - uniqParQ <- sort(unique(c(masterBlock))) - uniqParQ2 <- uniqParQ[which(uniqParQ > 0)] - concealnewQ <- (max(uniqParQ2) + 1):(max(uniqParQ2) + length(uniqParQ2)) - concealnewQMatr <- masterBlock - for (I in seq_along(uniqParQ2)) { - uniqParQ2 - concealnewQMatr[concealnewQMatr == uniqParQ2[I]] <- concealnewQ[I] + for (i in 1:ncol(Q)) { + for (j in 1:nrow(Q)) { + if (i == j) next + + obs_trait <- (c(i, j) - 1) %% n_obs_traits + 1 + conc_trait <- (c(i, j) - 1) %/% n_obs_traits + 1 + move_obs <- abs(diff(obs_trait)) + move_conc <- abs(diff(conc_trait)) + + new_rate <- 0 + + if (move_obs > 0 && move_conc == 0) { + + new_rate <- masterBlock[obs_trait[1], obs_trait[2]] + + } else if (move_conc > 0 && move_obs == 0) { + new_rate <- masterBlock[conc_trait[1], conc_trait[2]] + if (new_rate > 0) { + conc_index <- which(uniq_obs_rates == new_rate) + new_rate <- new_conceal_rates[conc_index] } - - Q <- create_q_matrix_int(masterBlock, - concealnewQMatr, - ntraits, - diff.conceal) - } - uniq_traits <- unique(traits) - uniq_traits <- uniq_traits[!is.na(uniq_traits)] - if (is.numeric(uniq_traits)) { - uniq_traits <- sort(uniq_traits) + } + + Q[i, j] <- new_rate } - - all_names <- get_state_names(state_names = uniq_traits, - num_concealed_states = length(uniq_traits)) - colnames(Q) <- all_names - rownames(Q) <- all_names - return(Q) + } + uniq_traits <- unique(traits) + uniq_traits <- uniq_traits[!is.na(uniq_traits)] + if (is.numeric(uniq_traits)) { + uniq_traits <- sort(uniq_traits) + } + + all_names <- get_state_names(state_names = uniq_traits, + num_concealed_states = num_concealed_states) + colnames(Q) <- all_names + rownames(Q) <- all_names + diag(Q) <- NA + return(Q) } + #' @title Data checking and trait sorting #' In preparation for likelihood calculation, it orders trait data according #' the tree tips @@ -427,8 +396,7 @@ check_traits <- function(traits, sampling_fraction) { if (all(sort(unique(as.vector(traits))) == sort(unique(traits[, 1]))) == FALSE) { - stop( - "Check your trait argument; if you have more than one column, + stop("Check your trait argument; if you have more than one column, make sure all your states are included in the first column." ) } @@ -460,9 +428,10 @@ check_root_state_weight <- function(root_state_weight, traits) { } else { if (any(root_state_weight == "maddison_weights" | root_state_weight == "equal_weights" | - root_state_weight == "proper_weights") == FALSE) { + root_state_weight == "proper_weights" | + root_state_weight == "stationary_weights") == FALSE) { stop("The root_state_weight must be any of - maddison_weights, equal_weights, or proper_weights.") + maddison_weights, equal_weights, proper_weights or stationary_weights.") } } } @@ -871,12 +840,16 @@ build_states <- function(phy, num_unique_traits = NULL, first_time = FALSE, traitStates = NULL) { + if (length(phy$tip.label) == 1) { + if (length(traits) > 1) { + traits <- matrix(traits, nrow = 1, ncol = length(traits)) + } + } if (!is.matrix(traits)) { traits <- matrix(traits, nrow = length(traits), ncol = 1, byrow = FALSE) } - if (length(phy$tip.label) != nrow(traits)) { - stop("Number of species in the tree must be the same as in the trait file") + stop("Number of species in the tree must be the same as in the trait file") } # if there are traits that are not in the observed tree, @@ -895,7 +868,7 @@ build_states <- function(phy, obs_traits <- unique(traits[, 1]) obs_traits <- obs_traits[!is.na(obs_traits)] if (sum(obs_traits %in% traitStates) != length(obs_traits)) { - stop("Tip traits are not in idparslist") + warning("Tip traits are not in idparslist") } nb_tip <- ape::Ntip(phy) @@ -939,6 +912,29 @@ build_initStates_time <- function(phy, num_unique_traits = NULL, first_time = FALSE, traitStates = NULL) { + + if (length(phy$tip.label) == 1) { + fake_phy <- ape::rphylo(n = 2, birth = 1, death = 0) + fake_phy$edge.length[1:2] <- phy$edge.length[1] + + fake_traits <- matrix(data = rep(traits, 2), + nrow = 2) + + states <- build_states(fake_phy, + fake_traits, + num_concealed_states, + sampling_fraction, + is_complete_tree, + mus, + num_unique_traits, + first_time, + traitStates) + phy$node.label <- NULL + states <- states[1, ] # only retain entry for one tip + forTime <- c(0, phy$edge.length) + ances <- NULL # this doesn't exist in a singleton tree + } else { + states <- build_states(phy, traits, num_concealed_states, @@ -948,12 +944,14 @@ build_initStates_time <- function(phy, num_unique_traits, first_time, traitStates) - phy$node.label <- NULL - split_times <- sort(event_times(phy), decreasing = FALSE) - ances <- as.numeric(names(split_times)) - - forTime <- cbind(phy$edge, phy$edge.length) - + + phy$node.label <- NULL + split_times <- sort(event_times(phy), decreasing = FALSE) + ances <- as.numeric(names(split_times)) + + forTime <- cbind(phy$edge, phy$edge.length) + } + return(list( states = states, ances = ances, @@ -968,7 +966,8 @@ get_weight_states <- function(root_state_weight, lambdas, nodeM, d, - is_cla = FALSE) { + is_cla = FALSE, + Q) { if (is.numeric(root_state_weight)) { weight_states <- rep(root_state_weight / num_concealed_states, @@ -998,8 +997,35 @@ get_weight_states <- function(root_state_weight, if (root_state_weight == "equal_weights") { weight_states <- rep(1 / length(mergeBranch), length(mergeBranch)) } + + if (root_state_weight == "stationary_weights") { + diag(Q) <- 0 + diag(Q) <- -rowSums(Q) + pi <- pracma::null(t(Q)) + diff <- 1 + dimpi2 <- dim(pi)[[2]] + if (dimpi2 > 1) { + warning('Null space of transition matrix is multidimensional; + the most even dimension is used for weighing the likelihood.') + } + for(i in 1:dimpi2) { + if (pi[which.max(abs(pi[,i])),i] < 0) { + pi[,i] <- -pi[,i] + } + if (any(pi[,i] < 0) && max(abs(pi[which(pi[,i] < 0),i])) > 1E-10) { + warning('Substantial negative weights detected.') + } + pi[which(pi[,i] < 0),i] <- 0 + pi[,i] <- pi[,i]/sum(pi[,i]) + diff_new <- abs(max(pi[,i])) - abs(min(pi[,i])) + if (diff_new < diff) { + diff <- diff_new + i_choice <- i + } + } + weight_states <- pi[,i_choice] + } } - return(weight_states) } @@ -1138,7 +1164,7 @@ check_ml_conditions <- function(traits, idparsfix, parsfix) { if (is.matrix(traits)) { - warning("you are setting a model where some species have more + warning("You are setting up a model where some species have more than one trait state") } @@ -1170,7 +1196,7 @@ check_ml_conditions <- function(traits, #' @keywords internal get_trait_states <- function(idparslist, num_concealed_states, - display_warning = TRUE) { + display_warning = FALSE) { trait_names <- names(idparslist[[1]]) if (is.null(trait_names)) trait_names <- names(idparslist[[2]]) if (is.null(trait_names)) trait_names <- colnames(idparslist[[3]]) @@ -1319,4 +1345,4 @@ get_root_state <- function(ancestral_states, phy, mus, d) { names(root_state) <- names(mus) return(root_state) -} \ No newline at end of file +} diff --git a/man/cla_secsse_loglik.Rd b/man/cla_secsse_loglik.Rd index 90a16c3..76d7dd5 100644 --- a/man/cla_secsse_loglik.Rd +++ b/man/cla_secsse_loglik.Rd @@ -50,10 +50,7 @@ to the number of examined states in the dataset.} \code{"proper_cond"} (default). For details, see vignette.} \item{root_state_weight}{the method to weigh the states: -\code{"maddison_weights"}, \code{"proper_weights"} (default) or \code{"equal_weights"}. -It can also be specified for the root state: the vector \code{c(1, 0, 0)} -indicates state 1 was the root state. When -using a \code{multiPhylo} object, root_state_weight should be list where each +\code{"maddison_weights"}, \code{"proper_weights"} (default), \verb{"equal_weights"'. or }"stationary_weights"\verb{It can also be specified for the root state: the vector}c(1, 0, 0)\verb{ indicates state 1 was the root state. When using a}multiPhylo` object, root_state_weight should be list where each entry in the list corresponds to the root_state_weight for each tree.} \item{sampling_fraction}{vector that states the sampling proportion per diff --git a/man/cla_secsse_ml.Rd b/man/cla_secsse_ml.Rd index 011441f..cf3af1d 100644 --- a/man/cla_secsse_ml.Rd +++ b/man/cla_secsse_ml.Rd @@ -63,10 +63,7 @@ to be estimated.} \code{"proper_cond"} (default). For details, see vignette.} \item{root_state_weight}{the method to weigh the states: -\code{"maddison_weights"}, \code{"proper_weights"} (default) or \code{"equal_weights"}. -It can also be specified for the root state: the vector \code{c(1, 0, 0)} -indicates state 1 was the root state. When -using a \code{multiPhylo} object, root_state_weight should be list where each +\code{"maddison_weights"}, \code{"proper_weights"} (default), \verb{"equal_weights"'. or }"stationary_weights"\verb{It can also be specified for the root state: the vector}c(1, 0, 0)\verb{ indicates state 1 was the root state. When using a}multiPhylo` object, root_state_weight should be list where each entry in the list corresponds to the root_state_weight for each tree.} \item{sampling_fraction}{vector that states the sampling proportion per diff --git a/man/cla_secsse_ml_func_def_pars.Rd b/man/cla_secsse_ml_func_def_pars.Rd index 9403ba0..fd3d577 100644 --- a/man/cla_secsse_ml_func_def_pars.Rd +++ b/man/cla_secsse_ml_func_def_pars.Rd @@ -81,10 +81,7 @@ example.} \code{"proper_cond"} (default). For details, see vignette.} \item{root_state_weight}{the method to weigh the states: -\code{"maddison_weights"}, \code{"proper_weights"} (default) or \code{"equal_weights"}. -It can also be specified for the root state: the vector \code{c(1, 0, 0)} -indicates state 1 was the root state. When -using a \code{multiPhylo} object, root_state_weight should be list where each +\code{"maddison_weights"}, \code{"proper_weights"} (default), \verb{"equal_weights"'. or }"stationary_weights"\verb{It can also be specified for the root state: the vector}c(1, 0, 0)\verb{ indicates state 1 was the root state. When using a}multiPhylo` object, root_state_weight should be list where each entry in the list corresponds to the root_state_weight for each tree.} \item{sampling_fraction}{vector that states the sampling proportion per diff --git a/man/default_params_doc.Rd b/man/default_params_doc.Rd index 30a5403..92daa28 100644 --- a/man/default_params_doc.Rd +++ b/man/default_params_doc.Rd @@ -109,10 +109,7 @@ fixed factors, so use a constant within \code{functions_defining_params}.} \code{"proper_cond"} (default). For details, see vignette.} \item{root_state_weight}{the method to weigh the states: -\code{"maddison_weights"}, \code{"proper_weights"} (default) or \code{"equal_weights"}. -It can also be specified for the root state: the vector \code{c(1, 0, 0)} -indicates state 1 was the root state. When -using a \code{multiPhylo} object, root_state_weight should be list where each +\code{"maddison_weights"}, \code{"proper_weights"} (default), \verb{"equal_weights"'. or }"stationary_weights"\verb{It can also be specified for the root state: the vector}c(1, 0, 0)\verb{ indicates state 1 was the root state. When using a}multiPhylo` object, root_state_weight should be list where each entry in the list corresponds to the root_state_weight for each tree.} \item{sampling_fraction}{vector that states the sampling proportion per diff --git a/man/plot_state_exact.Rd b/man/plot_state_exact.Rd index 3562d2a..157cdd7 100644 --- a/man/plot_state_exact.Rd +++ b/man/plot_state_exact.Rd @@ -48,10 +48,7 @@ entry in the list corresponds to the sampling proportion for each tree.} \code{"proper_cond"} (default). For details, see vignette.} \item{root_state_weight}{the method to weigh the states: -\code{"maddison_weights"}, \code{"proper_weights"} (default) or \code{"equal_weights"}. -It can also be specified for the root state: the vector \code{c(1, 0, 0)} -indicates state 1 was the root state. When -using a \code{multiPhylo} object, root_state_weight should be list where each +\code{"maddison_weights"}, \code{"proper_weights"} (default), \verb{"equal_weights"'. or }"stationary_weights"\verb{It can also be specified for the root state: the vector}c(1, 0, 0)\verb{ indicates state 1 was the root state. When using a}multiPhylo` object, root_state_weight should be list where each entry in the list corresponds to the root_state_weight for each tree.} \item{is_complete_tree}{logical specifying whether or not a tree with all its diff --git a/man/q_doubletrans.Rd b/man/q_doubletrans.Rd index 88cc5d9..903ab49 100644 --- a/man/q_doubletrans.Rd +++ b/man/q_doubletrans.Rd @@ -5,7 +5,7 @@ \title{Basic Qmatrix Sets a Q matrix where double transitions are not allowed} \usage{ -q_doubletrans(traits, masterBlock, diff.conceal) +q_doubletrans(traits, masterBlock, diff.conceal, num_concealed_states = NULL) } \arguments{ \item{traits}{vector with trait states for each tip in the phylogeny. The @@ -23,29 +23,31 @@ different. E.g. that the transition rates for the concealed states are different from the transition rates for the examined states. Normally it should be \code{FALSE} in order to avoid having a huge number of parameters.} + +\item{num_concealed_states}{number of concealed states, generally equivalent +to the number of examined states in the dataset.} } \value{ Q matrix that includes both examined and concealed states, it should be declared as the third element of idparslist. } \description{ -This function expands the Q_matrix, but it does so assuming -that the number of concealed traits is equal to the number of examined -traits, if you have a different number, you should consider looking at -the function \code{\link[=expand_q_matrix]{expand_q_matrix()}}. +This function expands the Q_matrix. If the number of concealed +states is not explicitly set by the user, it is assumed to be identical +to the number of observed states. } \examples{ -traits <- sample(c(0,1,2), 45,replace = TRUE) #get some traits +traits <- sample(c(0, 1, 2), 45, replace = TRUE) # get some traits # For a three-state trait -masterBlock <- matrix(99,ncol = 3,nrow = 3,byrow = TRUE) +masterBlock <- matrix(99, ncol = 3, nrow = 3, byrow = TRUE) diag(masterBlock) <- NA -masterBlock[1,2] <- 6 -masterBlock[1,3] <- 7 -masterBlock[2,1] <- 8 -masterBlock[2,3] <- 9 -masterBlock[3,1] <- 10 -masterBlock[3,2] <- 11 -myQ <- q_doubletrans(traits,masterBlock,diff.conceal = FALSE) +masterBlock[1, 2] <- 6 +masterBlock[1, 3] <- 7 +masterBlock[2, 1] <- 8 +masterBlock[2, 3] <- 9 +masterBlock[3, 1] <- 10 +masterBlock[3, 2] <- 11 +myQ <- q_doubletrans(traits, masterBlock, diff.conceal = FALSE) # now, it can replace the Q matrix from id_paramPos num_concealed_states <- 3 param_posit <- id_paramPos(traits,num_concealed_states) diff --git a/man/secsse_loglik.Rd b/man/secsse_loglik.Rd index fbab495..63774a8 100755 --- a/man/secsse_loglik.Rd +++ b/man/secsse_loglik.Rd @@ -50,10 +50,7 @@ to the number of examined states in the dataset.} \code{"proper_cond"} (default). For details, see vignette.} \item{root_state_weight}{the method to weigh the states: -\code{"maddison_weights"}, \code{"proper_weights"} (default) or \code{"equal_weights"}. -It can also be specified for the root state: the vector \code{c(1, 0, 0)} -indicates state 1 was the root state. When -using a \code{multiPhylo} object, root_state_weight should be list where each +\code{"maddison_weights"}, \code{"proper_weights"} (default), \verb{"equal_weights"'. or }"stationary_weights"\verb{It can also be specified for the root state: the vector}c(1, 0, 0)\verb{ indicates state 1 was the root state. When using a}multiPhylo` object, root_state_weight should be list where each entry in the list corresponds to the root_state_weight for each tree.} \item{sampling_fraction}{vector that states the sampling proportion per diff --git a/man/secsse_loglik_eval.Rd b/man/secsse_loglik_eval.Rd index ad6b933..f300c22 100644 --- a/man/secsse_loglik_eval.Rd +++ b/man/secsse_loglik_eval.Rd @@ -46,10 +46,7 @@ to the number of examined states in the dataset.} \code{"proper_cond"} (default). For details, see vignette.} \item{root_state_weight}{the method to weigh the states: -\code{"maddison_weights"}, \code{"proper_weights"} (default) or \code{"equal_weights"}. -It can also be specified for the root state: the vector \code{c(1, 0, 0)} -indicates state 1 was the root state. When -using a \code{multiPhylo} object, root_state_weight should be list where each +\code{"maddison_weights"}, \code{"proper_weights"} (default), \verb{"equal_weights"'. or }"stationary_weights"\verb{It can also be specified for the root state: the vector}c(1, 0, 0)\verb{ indicates state 1 was the root state. When using a}multiPhylo` object, root_state_weight should be list where each entry in the list corresponds to the root_state_weight for each tree.} \item{sampling_fraction}{vector that states the sampling proportion per diff --git a/man/secsse_ml.Rd b/man/secsse_ml.Rd index 38e82b7..ae20f31 100644 --- a/man/secsse_ml.Rd +++ b/man/secsse_ml.Rd @@ -63,10 +63,7 @@ to be estimated.} \code{"proper_cond"} (default). For details, see vignette.} \item{root_state_weight}{the method to weigh the states: -\code{"maddison_weights"}, \code{"proper_weights"} (default) or \code{"equal_weights"}. -It can also be specified for the root state: the vector \code{c(1, 0, 0)} -indicates state 1 was the root state. When -using a \code{multiPhylo} object, root_state_weight should be list where each +\code{"maddison_weights"}, \code{"proper_weights"} (default), \verb{"equal_weights"'. or }"stationary_weights"\verb{It can also be specified for the root state: the vector}c(1, 0, 0)\verb{ indicates state 1 was the root state. When using a}multiPhylo` object, root_state_weight should be list where each entry in the list corresponds to the root_state_weight for each tree.} \item{sampling_fraction}{vector that states the sampling proportion per diff --git a/man/secsse_ml_func_def_pars.Rd b/man/secsse_ml_func_def_pars.Rd index 20566e3..e63d5af 100644 --- a/man/secsse_ml_func_def_pars.Rd +++ b/man/secsse_ml_func_def_pars.Rd @@ -78,10 +78,7 @@ example.} \code{"proper_cond"} (default). For details, see vignette.} \item{root_state_weight}{the method to weigh the states: -\code{"maddison_weights"}, \code{"proper_weights"} (default) or \code{"equal_weights"}. -It can also be specified for the root state: the vector \code{c(1, 0, 0)} -indicates state 1 was the root state. When -using a \code{multiPhylo} object, root_state_weight should be list where each +\code{"maddison_weights"}, \code{"proper_weights"} (default), \verb{"equal_weights"'. or }"stationary_weights"\verb{It can also be specified for the root state: the vector}c(1, 0, 0)\verb{ indicates state 1 was the root state. When using a}multiPhylo` object, root_state_weight should be list where each entry in the list corresponds to the root_state_weight for each tree.} \item{sampling_fraction}{vector that states the sampling proportion per diff --git a/man/secsse_single_branch_loglik.Rd b/man/secsse_single_branch_loglik.Rd index 9ad536a..e397f5d 100644 --- a/man/secsse_single_branch_loglik.Rd +++ b/man/secsse_single_branch_loglik.Rd @@ -50,10 +50,7 @@ to the number of examined states in the dataset.} \code{"proper_cond"} (default). For details, see vignette.} \item{root_state_weight}{the method to weigh the states: -\code{"maddison_weights"}, \code{"proper_weights"} (default) or \code{"equal_weights"}. -It can also be specified for the root state: the vector \code{c(1, 0, 0)} -indicates state 1 was the root state. When -using a \code{multiPhylo} object, root_state_weight should be list where each +\code{"maddison_weights"}, \code{"proper_weights"} (default), \verb{"equal_weights"'. or }"stationary_weights"\verb{It can also be specified for the root state: the vector}c(1, 0, 0)\verb{ indicates state 1 was the root state. When using a}multiPhylo` object, root_state_weight should be list where each entry in the list corresponds to the root_state_weight for each tree.} \item{sampling_fraction}{vector that states the sampling proportion per diff --git a/src/Makevars.win b/src/Makevars.win index adf8485..11e7b47 100755 --- a/src/Makevars.win +++ b/src/Makevars.win @@ -1,4 +1,3 @@ CXX_STD = CXX17 -PKG_CXXFLAGS += -DRCPP_PARALLEL_USE_TBB=1 PKG_CPPFLAGS = -D_HAS_AUTO_PTR_ETC=0 PKG_LIBS += $(shell "${R_HOME}/bin${R_ARCH_BIN}/Rscript.exe" -e "RcppParallel::RcppParallelLibs()") \ No newline at end of file diff --git a/src/config.h b/src/config.h index f573168..f5f4dd3 100755 --- a/src/config.h +++ b/src/config.h @@ -8,6 +8,8 @@ #ifndef SRC_CONFIG_H_ #define SRC_CONFIG_H_ +#include "tbb_stub.h" + // Special case to make use of some steppers that would include // boost/get_pointer.hpp #ifndef BOOST_NO_AUTO_PTR diff --git a/src/mutable_dist.h b/src/mutable_dist.h index 02c7577..0d04367 100644 --- a/src/mutable_dist.h +++ b/src/mutable_dist.h @@ -9,8 +9,7 @@ #define mutable_dist_h #include - - +#include struct mutable_distribution { @@ -50,6 +49,4 @@ struct mutable_distribution { std::uniform_real_distribution<> unif_dist = std::uniform_real_distribution<>(0, 1.0); }; - - #endif /* mutable_dist_h */ diff --git a/src/secsse_loglik.cpp b/src/secsse_loglik.cpp index 32904ed..130716c 100755 --- a/src/secsse_loglik.cpp +++ b/src/secsse_loglik.cpp @@ -16,15 +16,6 @@ namespace secsse { - // probably the cleanest way to retrieve RcppParallel's concurrency setting - // set by RcppParallel::setThreadOptions(numThreads) - size_t get_rcpp_num_threads() { - auto* nt_env = std::getenv("RCPP_PARALLEL_NUM_THREADS"); - return (nullptr == nt_env) - ? tbb::task_arena::automatic // -1 - : static_cast(std::atoi(nt_env)); - } - template Rcpp::List calc_ll(std::unique_ptr od, const Rcpp::IntegerVector& ances, diff --git a/src/secsse_loglik.h b/src/secsse_loglik.h index 80cebdd..915675f 100755 --- a/src/secsse_loglik.h +++ b/src/secsse_loglik.h @@ -6,6 +6,7 @@ #pragma once +#include "config.h" #include #include #include @@ -16,11 +17,6 @@ namespace secsse { - // retreives value set by RcppParallel::setThreadOptions(numThreads) - // or tbb::task_arena::automatic if missing. - size_t get_rcpp_num_threads(); - - using state_ptr = std::vector*; // Models of 'integration_node` diff --git a/src/secsse_sim.cpp b/src/secsse_sim.cpp index caab4d5..ab4530d 100755 --- a/src/secsse_sim.cpp +++ b/src/secsse_sim.cpp @@ -123,6 +123,7 @@ Rcpp::List secsse_sim_cpp(const std::vector& m_R, } } if (sim.run_info == done) { + tracker[ done ]++; break; } else { tracker[ sim.run_info ]++; @@ -183,7 +184,7 @@ Rcpp::List secsse_sim_cpp(const std::vector& m_R, } catch (const char* msg) { Rcpp::Rcout << msg << std::endl; } catch(...) { - ::Rf_error("c++ exception (unknown reason)"); + Rcpp::stop("c++ exception (unknown reason)"); } return NA_REAL; } diff --git a/src/tbb_stub.h b/src/tbb_stub.h new file mode 100644 index 0000000..b68c622 --- /dev/null +++ b/src/tbb_stub.h @@ -0,0 +1,91 @@ +// Copyright (c) 2026, Hanno Hildenbrandt +// +// Distributed under the Boost Software License, Version 1.0. (See +// accompanying file LICENSE_1_0.txt or copy at +// http://www.boost.org/LICENSE_1_0.txt) + +// bare minimal tbb-stub header. +// just enough to to let *this* package pass the +// RCPP_PARALLEL_USE_TBB shenanigan on Alpine Linux + +#pragma once + +#include +#include +#include // pull RCPP_PARALLEL_USE_TBB + + +#if RCPP_PARALLEL_USE_TBB == 0 + +// everything looks so single-threaded here :( + +#include + +namespace tbb { + + namespace task_arena { + + constexpr size_t automatic = size_t(-1); + + } // namespace task_arena + + + class global_control { + public: + enum parameter { + max_allowed_parallelism, + thread_stack_size, + terminate_on_exception + }; + + global_control(parameter /*p*/, size_t /*value*/) {} + ~global_control() {}; + static size_t active_value(parameter /*param*/); // undefined + }; + + + template + inline void parallel_for_each( InputIterator first, InputIterator last, Body&& body ) { + std::for_each(first, last, std::forward(body)); + } + + + template + inline void parallel_for(Index first, Index last, const Func f) { + for (; first != last; ++first) { + f(first); + } + } + + + template + inline void parallel_for(Index first, Index last, Index step, const Func f) { + for (; first != last; first += step) { + f(first); + } + } + + +} // namespce tbb + + +// function name is lying. +inline size_t get_rcpp_num_threads() { + return 1; +} + + +#else // if RCPP_PARLLEL_USE_TBB = 0 + + +// probably the cleanest way to retrieve RcppParallel's concurrency setting +// set by RcppParallel::setThreadOptions(numThreads) +inline size_t get_rcpp_num_threads() { + auto* nt_env = std::getenv("RCPP_PARALLEL_NUM_THREADS"); + return (nullptr == nt_env) + ? tbb::task_arena::automatic // -1 + : static_cast(std::atoi(nt_env)); +} + + +#endif diff --git a/tests/testthat/test_lambda_setup.R b/tests/testthat/test_lambda_setup.R index 86baa66..67eed4f 100644 --- a/tests/testthat/test_lambda_setup.R +++ b/tests/testthat/test_lambda_setup.R @@ -136,7 +136,60 @@ test_that("test q_doubletrans", { a1 <- q_doubletrans(traits, masterBlock, diff.conceal = FALSE) a2 <- q_doubletrans(traits, masterBlock, diff.conceal = TRUE) - a1 <- unique(as.vector(a1)) - a2 <- unique(as.vector(a2)) - testthat::expect_gt(length(a2), length(a1)) + a1v <- unique(as.vector(a1)) + a2v <- unique(as.vector(a2)) + testthat::expect_gt(length(a2v), length(a1v)) + + # test deprecated function + testthat::expect_warning( + a1_1 <- secsse::expand_q_matrix(masterBlock, + num_concealed_states = num_concealed_states, + diff.conceal = FALSE), + "This function is deprecated, please use q_doubletrans, piping result + to q_doubletrans, this may introduce inaccuracies." + ) + # checking as vector, because expand_q_matrix mangles the trait names + testthat::expect_equal(as.vector(a1), as.vector(a1_1)) + + # trigger alternative route, I don't know why this route exists + #masterBlock <- matrix(5.1, ncol = 3, nrow = 3, byrow = TRUE) + #diag(masterBlock) <- NA + #a1 <- q_doubletrans(traits, masterBlock, diff.conceal = FALSE) + #a2 <- q_doubletrans(traits, masterBlock, diff.conceal = TRUE) + + #a1 <- unique(as.vector(a1)) + #a2 <- unique(as.vector(a2)) + #testthat::expect_gt(length(a2), length(a1)) + + traits <- c(2, 0, 1, 0, 2, 0, 1, 2, 2, 0) + num_concealed_states <- 3 + masterBlock <- matrix(5, ncol = 3, nrow = 3, byrow = TRUE) + diag(masterBlock) <- NA + a1 <- q_doubletrans(traits, masterBlock, diff.conceal = FALSE) + a2 <- q_doubletrans(traits, masterBlock, diff.conceal = TRUE) + testthat::expect_equal(max(a1, na.rm = TRUE), 5) + testthat::expect_equal(max(a2, na.rm = TRUE), 6) + + # less trivial example + traits <- c(1, 2, 3) + masterBlock <- matrix(c(0, 1, 0, + 2, 0, 3, + 0, 4, 0), ncol = 3, nrow = 3, byrow = TRUE) + diag(masterBlock) <- NA + a1 <- q_doubletrans(traits, masterBlock, diff.conceal = FALSE) + a2 <- q_doubletrans(traits, masterBlock, diff.conceal = TRUE) + testthat::expect_equal(max(a1, na.rm = TRUE), 4) + testthat::expect_equal(max(a2, na.rm = TRUE), 8) + + + a2 <- q_doubletrans(traits, masterBlock, diff.conceal = TRUE, + num_concealed_states = 2) + ref_mat <- matrix(c(0, 1, 0, 5, 0, 0, + 2, 0, 3, 0, 5 ,0, + 0, 4, 0, 0, 0, 5, + 6, 0, 0, 0, 1, 0, + 0, 6, 0, 2, 0, 3, + 0, 0, 6, 0, 4, 0), nrow = 6, byrow = TRUE) + diag(ref_mat) <- NA + testthat::expect_equal(as.vector(a2), as.vector(ref_mat)) }) \ No newline at end of file diff --git a/tests/testthat/test_multiphylo.R b/tests/testthat/test_multiphylo.R index 36f32a3..86e5a85 100644 --- a/tests/testthat/test_multiphylo.R +++ b/tests/testthat/test_multiphylo.R @@ -41,8 +41,8 @@ test_that("multi phylo", { cond = "no_cond", display_warning = FALSE) trees <- list() - trees[[1]]<- focal_tree - trees[[2]]<- focal_tree + trees[[1]] <- focal_tree + trees[[2]] <- focal_tree class(trees) <- "multiPhylo" @@ -71,6 +71,40 @@ test_that("multi phylo", { cond = "no_cond", display_warning = FALSE) testthat::expect_equal(2*res1, res3) + + # let's get root states + res4 <- secsse::cla_secsse_loglik(parameter = parslist, + phy = trees, + traits = trait_list, + num_concealed_states = num_concealed_states, + sampling_fraction = sf_list, + cond = "no_cond", + display_warning = FALSE, + return_root_state = TRUE) + testthat::expect_true(is.list(res4)) + testthat::expect_equal(length(res4$root_state), 2) + testthat::expect_equal(length(res4$root_state[[1]]), num_concealed_states ^ 2) + + # now, let's do multiphylo with a single branch tree + + fake_tree <- ape::rphylo(n = 2, birth = 1, death = 0) + fake_tree$edge.length <- c(1) + fake_tree$node.label <- NULL + fake_tree$edge <- as.matrix(fake_tree$edge[-2, ], nrow = 1) + fake_tree$tip.label <- fake_tree$tip.label[-2] + trees[[3]] <- fake_tree + trait_list[[3]] <- 1 + sf_list[[3]] <- c(1, 1) + + res5 <- secsse::cla_secsse_loglik(parameter = parslist, + phy = trees, + traits = trait_list, + num_concealed_states = num_concealed_states, + sampling_fraction = sf_list, + cond = "no_cond", + display_warning = FALSE) + # spec rate is 0.3, branch length is 1MY + testthat::expect_equal(res5 + params[1], res3, tol = 1e-2) }) test_that("multi phylo ML", { @@ -184,5 +218,95 @@ test_that("multi phylo ML", { testthat::expect_equal(3 * model_R$ML, multi_R$ML) testthat::expect_true(all.equal(model_R$MLpars, multi_R$MLpars)) + + # and now we add a single lineage tree: + fake_tree <- ape::rphylo(n = 2, birth = 1, death = 0) + fake_tree$edge.length <- c(1) + fake_tree$node.label <- NULL + fake_tree$edge <- as.matrix(fake_tree$edge[-2, ], nrow = 1) + fake_tree$tip.label <- fake_tree$tip.label[-2] + phylo_list[[4]] <- fake_tree + trait_list[[4]] <- 1 + sf_list[[4]] <- c(1, 1, 1) + + testthat::expect_message( + multi_R <- cla_secsse_ml( + phy = phylo_list, + traits = trait_list, + num_concealed_states = num_concealed_states, + idparslist = idparslist, + idparsopt = idparsopt, + initparsopt = initparsopt, + idparsfix = idparsfix, + parsfix = parsfix, + cond = cond, + root_state_weight = root_state_weight, + sampling_fraction = sf_list, + tol = tol, + maxiter = maxiter, + optimmethod = optimmethod, + num_cycles = 1, + verbose = FALSE) + ) + }) +test_that("multi phylo abuse", { + focal_tree <- ape::rphylo(n = 3, birth = 0.3 ,death = 0) + + traits <- c(1, 1, 1) + + num_concealed_states <- 2 + idparslist <- cla_id_paramPos(c(1, 2), num_concealed_states) + idparslist$lambdas[1, ] <- rep(1, 2) + idparslist[[2]][] <- 2 + masterBlock <- matrix(3, ncol = 2, nrow = 2, byrow = TRUE) + diag(masterBlock) <- NA + diff.conceal <- FALSE + idparslist[[3]] <- q_doubletrans(c(1, 2), masterBlock, diff.conceal) + idparslist[[1]] <- secsse::prepare_full_lambdas(c(1, 2), + num_concealed_states, + idparslist[[1]]) + + params <- c(0.3, 0.0, 0.0) # extinction and shifts to zero, to allow direct + # comparison + + lambdas <- secsse::fill_in(idparslist[[1]], params) + mus <- secsse::fill_in(idparslist[[2]], params) + q_mat <- secsse::fill_in(idparslist[[3]], params) + + parslist <- list() + parslist[[1]] <- lambdas + parslist[[2]] <- mus + parslist[[3]] <- q_mat + + sf <- c(1, 1) + + trees <- list() + trees[[1]] <- focal_tree + trees[[2]] <- focal_tree + + trait_list <- list() + trait_list[[1]] <- traits + trait_list[[2]] <- traits + + testthat::expect_error( + res2 <- secsse::cla_secsse_loglik(parameter = parslist, + phy = trees, + traits = trait_list, + num_concealed_states = num_concealed_states, + sampling_fraction = sf), + "when providing multiple phylogenies, make sure to use the multiPhylo class" + ) + + class(trees) <- "multiPhylo" + + testthat::expect_error( + res2 <- secsse::cla_secsse_loglik(parameter = parslist, + phy = trees, + traits = traits, + num_concealed_states = num_concealed_states, + sampling_fraction = sf), + "traits needs to be supplied as a list now that there are multiple phylogenies" + ) +}) \ No newline at end of file diff --git a/tests/testthat/test_plotting.R b/tests/testthat/test_plotting.R index 7b1f5de..f0aa55f 100644 --- a/tests/testthat/test_plotting.R +++ b/tests/testthat/test_plotting.R @@ -1,29 +1,28 @@ test_that("normal plotting", { - set.seed(5) - phy <- ape::rphylo(n = 4, birth = 1, death = 0) - traits <- c(0, 1, 1, 0) - params <- secsse::id_paramPos(c(0, 1), 2) - params[[1]][] <- c(0.2, 0.2, 0.1, 0.1) - params[[2]][] <- 0.01 - params[[3]][, ] <- 0.1 - diag(params[[3]]) <- NA - # Thus, we have for both, rates - # 0A, 1A, 0B and 1B. If we are interested in the posterior probability of - # trait 0,we have to provide a helper function that sums the probabilities of - # 0A and 0B, e.g.: - helper_function <- function(x) { - return(sum(x[c(5, 7)]) / sum(x)) # normalized by total sum, just in case. - } - testthat::expect_warning( - px <- plot_state_exact(parameters = params, - phy = phy, - traits = traits, - num_concealed_states = 2, - sampling_fraction = c(1, 1), - num_steps = 10, - prob_func = helper_function) - ) - testthat::expect_true(inherits(px, "ggplot")) + set.seed(5) + phy <- ape::rphylo(n = 4, birth = 1, death = 0) + traits <- c(0, 1, 1, 0) + params <- secsse::id_paramPos(c(0, 1), 2) + params[[1]][] <- c(0.2, 0.2, 0.1, 0.1) + params[[2]][] <- 0.01 + params[[3]][, ] <- 0.1 + diag(params[[3]]) <- NA + # Thus, we have for both, rates + # 0A, 1A, 0B and 1B. If we are interested in the posterior probability of + # trait 0,we have to provide a helper function that sums the probabilities of + # 0A and 0B, e.g.: + helper_function <- function(x) { + return(sum(x[c(5, 7)]) / sum(x)) # normalized by total sum, just in case. + } + px <- plot_state_exact(parameters = params, + phy = phy, + traits = traits, + num_concealed_states = 2, + sampling_fraction = c(1, 1), + num_steps = 10, + prob_func = helper_function) + + testthat::expect_true(inherits(px, "ggplot")) }) test_that("cla plotting", { @@ -38,9 +37,9 @@ test_that("cla plotting", { masterBlock <- matrix(5, ncol = 3, nrow = 3, byrow = TRUE) diag(masterBlock) <- NA diff.conceal <- FALSE - + idparslist[[3]] <- q_doubletrans(traits, masterBlock, diff.conceal) - + testthat::expect_output( startingpoint <- DDD::bd_ML(brts = ape::branching.times(phylotree)) ) @@ -55,7 +54,7 @@ test_that("cla plotting", { cond <- "proper_cond" root_state_weight <- "proper_weights" sampling_fraction <- c(1, 1, 1) - + model_R <- secsse::cla_secsse_ml( phy = phylotree, traits = traits, @@ -76,18 +75,17 @@ test_that("cla plotting", { helper_function <- function(x) { return(sum(x[c(10, 13, 16)]) / sum(x)) } - - testthat::expect_warning( - px <- secsse::plot_state_exact(parameters = model_R$MLpars, - phy = phylotree, - traits = traits, - num_concealed_states = + + px <- secsse::plot_state_exact(parameters = model_R$MLpars, + phy = phylotree, + traits = traits, + num_concealed_states = num_concealed_states, - sampling_fraction = sampling_fraction, - cond = cond, - root_state_weight = root_state_weight, - prob_func = helper_function) - ) - + sampling_fraction = sampling_fraction, + cond = cond, + root_state_weight = root_state_weight, + prob_func = helper_function) + + testthat::expect_true(inherits(px, "ggplot")) }) diff --git a/tests/testthat/test_root_state.R b/tests/testthat/test_root_state.R index 896d820..cb32c1a 100644 --- a/tests/testthat/test_root_state.R +++ b/tests/testthat/test_root_state.R @@ -40,7 +40,8 @@ test_that("root state cla LL", { qs = params[[3]], crown_age = 5, num_concealed_states = num_concealed_states, - sampling_fraction = c(1, 1)) + sampling_fraction = c(1, 1), + seed = 5) diff --git a/tests/testthat/test_secsse_sim.R b/tests/testthat/test_secsse_sim.R index 285c9c1..14f6e40 100644 --- a/tests/testthat/test_secsse_sim.R +++ b/tests/testthat/test_secsse_sim.R @@ -124,6 +124,20 @@ test_that("test secsse_sim 2", { vx <- geiger::is.extinct(focal_tree$phy) testthat::expect_true(length(vx) > 0) } + + # test verbose output + testthat::expect_output( + focal_tree2 <- secsse::secsse_sim(lambdas = lambda_p, + mus = mu_p, + qs = q_mat_p, + crown_age = 10, + num_concealed_states = 2, + max_spec = 100, + seed = 21, + verbose = TRUE, + drop_extinct = FALSE) + ) + testthat::expect_equal(focal_tree, focal_tree2) }) @@ -805,4 +819,151 @@ test_that("test secsse_sim sampling_fraction", { tolerance = 0.5) testthat::expect_equal(median(num_ones / num_zeros), sf, tolerance = 0.1) } + + + # now test something where due to the sampling fraction, it goes extinct + # This will actually generate multiple warnings, but I can not make them explicit + # all here. + testthat::expect_warning( + focal_tree <- secsse::secsse_sim(lambdas = lambda_p, + mus = mu_p, + qs = q_mat_p, + crown_age = crown_age, + num_concealed_states = 2, + sampling_fraction = c(0.0, 0.0), + seed = 42, + max_tries = 2) + ) }) + + +test_that("test secsse_sim root edge", { + + states <- c("S", "G") + + spec_matrix <- c("S", "S", "S", 1) + spec_matrix <- rbind(spec_matrix, c("G", "G", "G", 2)) + lambda_list <- secsse::create_lambda_list(state_names = states, + num_concealed_states = 2, + transition_matrix = spec_matrix, + model = "ETD") + + mu_vector <- secsse::create_mu_vector(state_names = states, + num_concealed_states = 2, + model = "ETD", + lambda_list = lambda_list) + + shift_matrix <- c("S", "G", 5) + shift_matrix <- rbind(shift_matrix, c("G", "S", 6)) + + q_matrix <- secsse::create_q_matrix(state_names = states, + num_concealed_states = 2, + shift_matrix = shift_matrix, + diff.conceal = TRUE) + idparslist <- list() + idparslist[[1]] <- lambda_list + idparslist[[2]] <- mu_vector + idparslist[[3]] <- q_matrix + + spec_S <- 0.1 + spec_G <- 0 + ext_S <- ext_G <- 0.0 + q_SG <- 0.0 + q_GS <- 0 + used_params <- c(spec_S, spec_G, ext_S, ext_G, q_SG, q_GS, 0, 0) + + crown_age_used <- 10 + + sim_lambda_list_etd <- secsse::fill_in(idparslist[[1]], used_params) + sim_mu_vector_etd <- secsse::fill_in(idparslist[[2]], used_params) + sim_q_matrix_etd <- secsse::fill_in(idparslist[[3]], used_params) + + + sim_tree <- secsse::secsse_sim(lambdas = sim_lambda_list_etd, + mus = sim_mu_vector_etd, + qs = sim_q_matrix_etd, + crown_age = crown_age_used, + num_concealed_states = 2, + conditioning = "none", + min_spec = 1, + init_state_probs = c("SA", "SB"), + start_at_crown = FALSE, + seed = 1) + testthat::expect_equal(unique(sim_tree$obs_traits), "S") + testthat::expect_true(!is.null(sim_tree$phy$root.edge)) + + # and now including extinction + + spec_S <- 0.2 + spec_G <- 0 + ext_S <- ext_G <- 0.15 + q_SG <- 0 + q_GS <- 0 + used_params <- c(spec_S, spec_G, ext_S, ext_G, q_SG, q_GS, 0, 0) + + crown_age_used <- 10 + + sim_lambda_list_etd <- secsse::fill_in(idparslist[[1]], used_params) + sim_mu_vector_etd <- secsse::fill_in(idparslist[[2]], used_params) + sim_q_matrix_etd <- secsse::fill_in(idparslist[[3]], used_params) + + + sim_tree2 <- secsse::secsse_sim(lambdas = sim_lambda_list_etd, + mus = sim_mu_vector_etd, + qs = sim_q_matrix_etd, + crown_age = crown_age_used, + num_concealed_states = 2, + conditioning = "none", + min_spec = 1, + init_state_probs = c("SA", "SB"), + start_at_crown = FALSE, + seed = 3) + testthat::expect_true(!is.null(sim_tree2$phy$root.edge)) + testthat::expect_equal(unique(sim_tree2$obs_traits), "S") + testthat::expect_gt(sim_tree2$extinct, sim_tree$extinct) + + # single lineage tree: + spec_S <- 0.01 + used_params <- c(spec_S, spec_G, ext_S, ext_G, q_SG, q_GS, 0, 0) + + sim_lambda_list_etd <- secsse::fill_in(idparslist[[1]], used_params) + + sim_tree <- secsse::secsse_sim(lambdas = sim_lambda_list_etd, + mus = sim_mu_vector_etd, + qs = sim_q_matrix_etd, + crown_age = crown_age_used, + num_concealed_states = 2, + conditioning = "none", + min_spec = 1, + init_state_probs = c("SA", "SB"), + start_at_crown = FALSE, + seed = 1) + testthat::expect_equal(length(sim_tree$phy$tip.label), 1) + testthat::expect_equal(sim_tree$phy$edge.length[1], crown_age_used) + + # force failure of simulation of tree + # single lineage tree: + ext_S <- 10 + used_params <- c(spec_S, spec_G, ext_S, ext_G, q_SG, q_GS, 0, 0) + + sim_mu_vector_etd <- secsse::fill_in(idparslist[[2]], used_params) + + testthat::expect_warning( + sim_tree <- secsse::secsse_sim(lambdas = sim_lambda_list_etd, + mus = sim_mu_vector_etd, + qs = sim_q_matrix_etd, + crown_age = crown_age_used, + num_concealed_states = 2, + conditioning = "none", + min_spec = 1, + init_state_probs = c("SA", "SB"), + start_at_crown = FALSE, + seed = 1, + max_tries = 100), + "Couldn't simulate a tree in enough tries, + try increasing max_tries" + ) + + testthat::expect_equal(sim_tree$status, "not enough tries") + testthat::expect_equal(sim_tree$phy, "ds") +}) \ No newline at end of file diff --git a/tests/testthat/test_weights.R b/tests/testthat/test_weights.R new file mode 100644 index 0000000..20734ed --- /dev/null +++ b/tests/testthat/test_weights.R @@ -0,0 +1,112 @@ +test_that("stationary weights are computed correctly", { + Sys.unsetenv("R_TESTS") + Q <- matrix(c(0,3,2,0), 2, 2, byrow = T) + weight_states <- get_weight_states(root_state_weight = 'stationary_weights', + num_concealed_states = NA, + mergeBranch = NA, + lambdas = NA, + nodeM = NA, + d = NA, + is_cla = NA, + Q = Q) + testthat::expect_equal(as.numeric(weight_states),c(0.4,0.6)) + + q <- matrix(0, ncol = 3, nrow = 3, byrow = TRUE) + q[2, 1] <- 1.4 + q[3, 1] <- 1.3 + q[1, 2] <- 0.7 + q[1, 3] <- 0.7 + + weight_states <- get_weight_states(root_state_weight = 'stationary_weights', + num_concealed_states = NA, + mergeBranch = NA, + lambdas = NA, + nodeM = NA, + d = NA, + is_cla = NA, + Q = q) + testthat::expect_equal(as.numeric(weight_states),c(0.4905660,0.2452830,0.2641509), tol = 1E-6) +}) + +test_that("secsse works with stationary weights", { + Sys.unsetenv("R_TESTS") + + utils::data("example_phy_GeoSSE", package = "secsse") + traits <- as.numeric(example_phy_GeoSSE$tip.state) + + lambdas <- list() + lambdas[[1]] <- matrix(0, ncol = 3, nrow = 3, byrow = TRUE) + lambdas[[2]] <- lambdas[[1]] + lambdas[[3]] <- lambdas[[1]] + lambdas[[1]][1, 1] <- 1.5 + lambdas[[2]][2, 2] <- 0.5 + lambdas[[3]][3, 3] <- 1 + + mus <- c(0.7, 0.7, 0.7) + + q <- matrix(0, ncol = 3, nrow = 3, byrow = TRUE) + q[2, 1] <- 1.4 + q[3, 1] <- 1.3 + q[1, 2] <- 0.7 + q[1, 3] <- 0.7 + + parameter <- list() + parameter[[1]] <- lambdas + parameter[[2]] <- mus + parameter[[3]] <- q + + num_concealed_states <- 3 + + num_modeled_traits <- ncol(q) / floor(num_concealed_states) + + setting_calculation <- build_initStates_time(phy = example_phy_GeoSSE, + traits = traits, + num_concealed_states = num_concealed_states, + sampling_fraction = c(1, 1, 1), + is_complete_tree = FALSE, + mus = mus, + num_unique_traits = num_modeled_traits, + first_time = TRUE) + states <- setting_calculation$states + d <- ncol(states) / 3 + new_states <- states[, c(1, 2, 3, 10, 11, 12, 19, 20, 21)] + states <- new_states + setting_calculation$states <- states + + cla_secsse_LL <- secsse::cla_secsse_loglik(parameter = parameter, + phy = example_phy_GeoSSE, + traits = traits, + num_concealed_states = 3, + cond = "proper_cond", + root_state_weight = "stationary_weights", + sampling_fraction = c(1, 1, 1), + setting_calculation = setting_calculation, + see_ancestral_states = FALSE, + loglik_penalty = 0, + is_complete_tree = FALSE, + num_threads = 1, + method = "odeint::runge_kutta_cash_karp54", + atol = 1e-8, + rtol = 1e-7) + + pars <- parameter + pars[[1]] <- c(lambdas[[1]][1, 1], lambdas[[2]][2, 2], lambdas[[3]][3, 3]) + + secsse_LL <- secsse::secsse_loglik(parameter = pars, + phy = example_phy_GeoSSE, + traits = traits, + num_concealed_states = num_concealed_states, + cond = "proper_cond", + root_state_weight = "stationary_weights", + sampling_fraction = c(1,1,1), + setting_calculation = setting_calculation, + see_ancestral_states = FALSE, + loglik_penalty = 0, + is_complete_tree = FALSE, + num_threads = 1, + atol = 1e-8, + rtol = 1e-7, + method = "odeint::runge_kutta_cash_karp54") + + testthat::expect_equal(secsse_LL, cla_secsse_LL,tolerance = 1e-5) +}) diff --git a/vignettes/plotting_states.Rmd b/vignettes/plotting_states.Rmd index 149c89d..28d0bfd 100644 --- a/vignettes/plotting_states.Rmd +++ b/vignettes/plotting_states.Rmd @@ -63,18 +63,19 @@ Here, the first four rows indicate the tip states, whilst the later three rows indicate the states at the internal nodes (with the last row indicating the root, in this case). The columns indicate the four extinction and four speciation rates, following the order in params[[1]] and params[[2]]. Thus, we have for both, rates -0A, 1A, 0B and 1B. If we are interested in the posterior probability of trait 0, +0A, 1A, 0B and 1B. Furthermore, the last four columns indicate the survival rates, +which are equal to 1 - extinction. If we are interested in the posterior probability of trait 0, we have to provide a helper function that sums the probabilities of 0A and 0B, e.g.: ```{r helper function} helper_function <- function(x) { - return(sum(x[c(5, 7)]) / sum(x)) # normalized by total sum, just in case. + return(sum(x[c(5, 7)]) / sum(x[1:8])) # normalized by total sum, just in case. } ``` We can now use this to plot this probability across the tree. There are two options for plotting: using the evaluations along the branches as used by the integration method, or evaluating the branch values at a specific number of intervals. -Using the explicit evaluations is more precies, but might be memory heavy. Usually, +Using the explicit evaluations is more precise, but might be memory heavy. Usually, using 10-100 evaluations per branch provides a very accurate approximation: ```{r exact} diff --git a/vignettes/starting_secsse.Rmd b/vignettes/starting_secsse.Rmd index 9b18315..dd4ecbc 100644 --- a/vignettes/starting_secsse.Rmd +++ b/vignettes/starting_secsse.Rmd @@ -97,23 +97,24 @@ parameters and constraints. #### Note on assigning ambiguity to taxon trait states -If the user wishes to assign a taxon to multiple trait states, because he/she is -unsure which state best describes the taxon, he/she can use `NA`. `NA` is used +If the user wishes to assign a taxon to multiple trait states, because she/he is +unsure which state best describes the taxon, she/he can use `NA`. `NA` is used when there is no information on possible state at all; for example when a state was not measured or a taxon is unavailable for inspection. `NA` means a taxon is equally likely to pertain to any state. In case the user does have some information, for example if a taxon can pertain to multiple states, or if there is uncertainty regarding state but one or multiple states can with certainty be excluded, secsse offers flexibility to handle ambiguity. In this case, the user -only needs to supply a trait file, with at least four columns, one for the taxon -name, and three for trait state. Below, we show an example of what the trait +only needs to supply a trait matrix, with at least N + 1 columns, one column for +the taxon name, and N equal to the number of unique states (3 in the example) . +Below, we show an example of what the trait info should be like (the column with species' names has been removed). If a taxon may pertain to trait state 1 or 3, but not to 2, the three columns should -have at least the values 1 and a 3, but never 2 (species in the third row). On +have at least the values 1 and a 3, but never 2 (species in the fourth row). On the other hand, the species in the fifth row can pertain to all states: the first column would have a 1, the second a 2, the third a 3 (although if you only have this type of ambiguity, it is easier to assign `NA` and use a single-column -data file). +setup). ```{r} # traits traits traits @@ -128,8 +129,7 @@ data file). ## Setting up an analysis To perform a Maximum Likelihood analysis, secsse makes use of the -function `DDD::optimize()`, which in turn, typically, uses the subplex -package to perform the Maximum Likelihood optimization. In such an +function `DDD::optimize()`. In such an analysis, we need to specify which parameters we want to optimize, which parameters to keep fix, and the initial values per parameter. We do so by providing the structure of the input parameters (e.g. in vector,