diff --git a/DESCRIPTION b/DESCRIPTION index c055a0e..7d7df1c 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.1 +Version: 3.7.0 License: GPL (>= 3) | file LICENSE Authors@R: c( person(given = "Leonel", diff --git a/NEWS.md b/NEWS.md index 2901193..8b928b0 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,8 @@ +# 3.7.0 +- Fixed a major bug when doing Maximum Likelihood including data sets that have +singleton trees (e.g. trees with only one tip) has been fixed. All analyses on +such a data set being done with a secsse version prior to this are incorrect. + # 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`. diff --git a/R/seccse_plot.R b/R/seccse_plot.R index 8666aec..2197648 100644 --- a/R/seccse_plot.R +++ b/R/seccse_plot.R @@ -141,6 +141,7 @@ plot_state_exact <- function(parameters, root_state_weight = "proper_weights", is_complete_tree = FALSE, method = "odeint::runge_kutta_cash_karp54", + num_threads = 1, atol = 1e-16, rtol = 1e-16, num_steps = 100, @@ -162,7 +163,8 @@ plot_state_exact <- function(parameters, is_complete_tree = is_complete_tree, atol = atol, rtol = rtol, - method = method) + method = method, + num_threads = num_threads) if (verbose) message("\nconverting collected likelihoods to graph positions:\n") diff --git a/man/plot_state_exact.Rd b/man/plot_state_exact.Rd index 157cdd7..6fbe0d3 100644 --- a/man/plot_state_exact.Rd +++ b/man/plot_state_exact.Rd @@ -14,6 +14,7 @@ plot_state_exact( root_state_weight = "proper_weights", is_complete_tree = FALSE, method = "odeint::runge_kutta_cash_karp54", + num_threads = 1, atol = 1e-16, rtol = 1e-16, num_steps = 100, @@ -60,6 +61,8 @@ extinct species is provided. If set to \code{TRUE}, it also assumes that all \code{"odeint::runge_kutta_dopri5"}, \code{"odeint::bulirsch_stoer"} and \code{"odeint::runge_kutta4"}. Default method is: \code{"odeint::runge_kutta_cash_karp54"}.} +\item{num_threads}{number of threads to be used. Default is one thread.} + \item{atol}{A numeric specifying the absolute tolerance of integration.} \item{rtol}{A numeric specifying the relative tolerance of integration.} diff --git a/src/config.h b/src/config.h index f5f4dd3..692789f 100755 --- a/src/config.h +++ b/src/config.h @@ -36,7 +36,7 @@ #define SECSSE_DEFAULT_EVAL_DTF 0.1 // Uncomment to enable nested parallelism. -// This feature may improve or may deterioate performance. +// This feature may improve or may deteriorate performance. // Default is disabled. //#define SECSSE_NESTED_PARALLELISM diff --git a/src/secsse_eval.cpp b/src/secsse_eval.cpp index ce48aaf..a0e9e23 100755 --- a/src/secsse_eval.cpp +++ b/src/secsse_eval.cpp @@ -38,14 +38,15 @@ Rcpp::List eval(std::unique_ptr od, tstates.emplace_back(states.row(i).begin(), states.row(i).end()); } const auto phy_edge = make_phy_edge_vector(rmatrix(forTime)); - auto inodes = find_inte_nodes(phy_edge, rvector(ances), tstates); + auto inodes = find_inte_nodes(phy_edge, rvector(ances), tstates, num_threads); auto integrator = Integrator(std::move(od), method, atol, rtol); - calc_ll(integrator, inodes, tstates); + calc_ll(integrator, inodes, tstates, num_threads); // integrate over each edge auto snodes = inodes_t(std::begin(inodes), std::end(inodes)); - tbb::parallel_for_each(std::begin(snodes), std::end(snodes), + tbb::task_arena(num_threads).execute([&] { + tbb::parallel_for_each(std::begin(snodes), std::end(snodes), [&](auto& snode) { #ifdef SECSSE_NESTED_PARALLELISM tbb::parallel_for(0, 2, [&](size_t i) { @@ -56,6 +57,7 @@ Rcpp::List eval(std::unique_ptr od, integrator(snode.desc[1], num_steps); #endif }); + }); // convert to Thijs's data layout: // rows of [ances, focal, t, [probs]] const size_t nrow = 2 * snodes.size() * (num_steps + 1); diff --git a/src/secsse_loglik.cpp b/src/secsse_loglik.cpp index 130716c..b8211c3 100755 --- a/src/secsse_loglik.cpp +++ b/src/secsse_loglik.cpp @@ -36,13 +36,13 @@ namespace secsse { tstates.emplace_back(states.row(i).begin(), states.row(i).end()); } const auto phy_edge = make_phy_edge_vector(rmatrix(forTime)); - auto inodes = find_inte_nodes(phy_edge, rvector(ances), tstates); + auto inodes = find_inte_nodes(phy_edge, rvector(ances), tstates, num_threads); calc_ll_res ll_res; if (use_normalization) { - ll_res = calc_ll(Integrator( std::move(od), method, atol, rtol), inodes, tstates); + ll_res = calc_ll(Integrator( std::move(od), method, atol, rtol), inodes, tstates, num_threads); } else { - ll_res = calc_ll(Integrator(std::move(od), method, atol, rtol), inodes, tstates); + ll_res = calc_ll(Integrator(std::move(od), method, atol, rtol), inodes, tstates, num_threads); } diff --git a/src/secsse_loglik.h b/src/secsse_loglik.h index 915675f..c279308 100755 --- a/src/secsse_loglik.h +++ b/src/secsse_loglik.h @@ -101,28 +101,29 @@ namespace secsse { } - inline inodes_t find_inte_nodes(const std::vector& phy_edge, rvector ances, std::vector>& states) { + inline inodes_t find_inte_nodes(const std::vector& phy_edge, rvector ances, std::vector>& states, const int num_threads) { auto res = inodes_t{ances.size()}; auto comp = [](auto& edge, size_t val) { return edge.n < val; }; - tbb::parallel_for(0, ances.size(), 1, [&](int i) { - const auto focal = ances[i]; - auto& inode = res[i]; - inode.state = &states[focal - 1]; - inode.state->clear(); // 'dirty' condition - auto it0 = std::lower_bound(std::begin(phy_edge), std::end(phy_edge), focal, comp); - auto it1 = std::lower_bound(it0 + 1, std::end(phy_edge), focal, comp); - // the next thingy is easy to overlook: the sequence matters for creating - // the 'merged' branch. imposes some pre-condition that is nowere to find :( - if (it0->m > it1->m) { - std::swap(it0, it1); - } - inode.desc[0] = { &states[it0->m - 1], it0->time }; - inode.desc[1] = { &states[it1->m - 1], it1->time }; + tbb::task_arena(num_threads).execute([&] { + tbb::parallel_for(0, ances.size(), 1, [&](int i) { + const auto focal = ances[i]; + auto& inode = res[i]; + inode.state = &states[focal - 1]; + inode.state->clear(); // 'dirty' condition + auto it0 = std::lower_bound(std::begin(phy_edge), std::end(phy_edge), focal, comp); + auto it1 = std::lower_bound(it0 + 1, std::end(phy_edge), focal, comp); + // the next thingy is easy to overlook: the sequence matters for creating + // the 'merged' branch. imposes some pre-condition that is nowhere to find :( + if (it0->m > it1->m) { + std::swap(it0, it1); + } + inode.desc[0] = { &states[it0->m - 1], it0->time }; + inode.desc[1] = { &states[it1->m - 1], it1->time }; + }); }); return res; } - template inline double normalize_loglik(RaIt first, RaIt last) { const auto sabs = std::accumulate(first, last, 0.0, [](const auto& s, const auto& x) { @@ -230,7 +231,8 @@ namespace secsse { template inline calc_ll_res calc_ll(const INTEGRATOR& integrator, inodes_t& inodes, - std::vector>& /* in/out */ states) + std::vector>& /* in/out */ states, + int num_threads) { const auto d = integrator.size(); auto is_dirty = [](const auto& inode) { @@ -238,8 +240,10 @@ namespace secsse { }; for (auto first = std::begin(inodes); first != std::end(inodes) ;) { auto last = std::partition(first, std::end(inodes), std::not_fn(is_dirty)); - tbb::parallel_for_each(first, last, [&](auto& inode) { - integrator(inode); + tbb::task_arena(num_threads).execute([&] { + tbb::parallel_for_each(first, last, [&](auto& inode) { + integrator(inode); + }); }); first = last; } diff --git a/src/tbb_stub.h b/src/tbb_stub.h index b68c622..006569a 100644 --- a/src/tbb_stub.h +++ b/src/tbb_stub.h @@ -65,7 +65,10 @@ namespace tbb { } } - + class task_arena { + task_arena(size_t /*num_threads*/) {}); + template void execute(const Func& f) { f(); } + }; } // namespce tbb diff --git a/vignettes/complete_tree.Rmd b/vignettes/complete_tree.Rmd index eaeb09f..1862464 100644 --- a/vignettes/complete_tree.Rmd +++ b/vignettes/complete_tree.Rmd @@ -136,7 +136,8 @@ complete_tree_ml_CR <- secsse_ml(phy = sim_tree_complete$phy, idparsfix = idparsfix, parsfix = initparsfix, sampling_fraction = sampling_fraction, - verbose = FALSE) + verbose = FALSE, + num_threads = 1) ``` Now we can see what our results look like. @@ -195,7 +196,8 @@ reconstructed_tree_ml <- secsse_ml(phy = sim_tree_reconstructed$phy, parsfix = initparsfix, sampling_fraction = sampling_fraction, verbose = FALSE, - is_complete_tree = FALSE) + is_complete_tree = FALSE, + num_threads = 1) ``` @@ -243,7 +245,7 @@ knitr::kable( ) ``` -We see that including extinct species results in a better esimation +We see that including extinct species results in a better estimation particularly of the extinction rate. This effect is especially noticeable if there are many extinct species present in the tree. Additionally, we see that the estimation of the transition rate from state 1 to diff --git a/vignettes/plotting_states.Rmd b/vignettes/plotting_states.Rmd index 28d0bfd..89f1b43 100644 --- a/vignettes/plotting_states.Rmd +++ b/vignettes/plotting_states.Rmd @@ -42,11 +42,12 @@ diag(params[[3]]) <- NA ll <- secsse::secsse_loglik(parameter = params, - phy = phy, - traits = traits, - num_concealed_states = 2, - see_ancestral_states = TRUE, - sampling_fraction = c(1, 1)) + phy = phy, + traits = traits, + num_concealed_states = 2, + see_ancestral_states = TRUE, + sampling_fraction = c(1, 1), + num_threads = 1) ll ``` @@ -84,7 +85,8 @@ secsse::plot_state_exact(parameters = params, traits = traits, num_concealed_states = 2, sampling_fraction = c(1, 1), - prob_func = helper_function) + prob_func = helper_function, + num_threads = 1) secsse::plot_state_exact(parameters = params, phy = phy, @@ -92,7 +94,8 @@ secsse::plot_state_exact(parameters = params, num_concealed_states = 2, sampling_fraction = c(1, 1), num_steps = 10, - prob_func = helper_function) + prob_func = helper_function, + num_threads = 1) secsse::plot_state_exact(parameters = params, phy = phy, @@ -100,7 +103,8 @@ secsse::plot_state_exact(parameters = params, num_concealed_states = 2, sampling_fraction = c(1, 1), num_steps = 100, - prob_func = helper_function) + prob_func = helper_function, + num_threads = 1) ``` # Using CLA secsse @@ -160,6 +164,7 @@ secsse::plot_state_exact(parameters = parameter, root_state_weight = "maddison_weights", is_complete_tree = FALSE, prob_func = helper_function, - num_steps = 10) + num_steps = 10, + num_threads = 1) ``` diff --git a/vignettes/starting_secsse.Rmd b/vignettes/starting_secsse.Rmd index dd4ecbc..13f2377 100644 --- a/vignettes/starting_secsse.Rmd +++ b/vignettes/starting_secsse.Rmd @@ -301,7 +301,8 @@ answ <- secsse::cla_secsse_ml(phy = phylo_vignette, idparsfix = idparsfix, parsfix = initparsfix, sampling_fraction = sampling_fraction, - verbose = FALSE) + verbose = FALSE, + num_threads = 1) ``` We can now extract several pieces of information from the returned @@ -419,7 +420,8 @@ answ <- secsse::cla_secsse_ml(phy = phylo_vignette, idparsfix = idparsfix, parsfix = initparsfix, sampling_fraction = sampling_fraction, - verbose = FALSE) + verbose = FALSE, + num_threads = 1) ML_CTD <- answ$ML CTD_par <- secsse::extract_par_vals(idparslist, answ$MLpars) ML_CTD @@ -522,7 +524,8 @@ answ <- secsse::cla_secsse_ml(phy = phylo_vignette, idparsfix = idparsfix, parsfix = initparsfix, sampling_fraction = sampling_fraction, - verbose = FALSE) + verbose = FALSE, + num_threads = 1) ML_CR <- answ$ML CR_par <- secsse::extract_par_vals(idparslist, answ$MLpars) ML_CR