Develop

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

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`.

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

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
 

src/secsse_eval.cpp

@@ -38,14 +38,15 @@ Rcpp::List eval(std::unique_ptr<ODE> od,
     tstates.emplace_back(states.row(i).begin(), states.row(i).end());
   }
   const auto phy_edge = make_phy_edge_vector(rmatrix<const double>(forTime));
-  auto inodes = find_inte_nodes(phy_edge, rvector<const int>(ances), tstates);
+  auto inodes = find_inte_nodes(phy_edge, rvector<const int>(ances), tstates, num_threads);
   auto integrator = Integrator<ODE, NORMALIZER>(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<storing::inode_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<ODE> 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);

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<const double>(forTime));
-    auto inodes = find_inte_nodes(phy_edge, rvector<const int>(ances), tstates);
+    auto inodes = find_inte_nodes(phy_edge, rvector<const int>(ances), tstates, num_threads);
 
     calc_ll_res ll_res;
     if (use_normalization) {
-      ll_res  = calc_ll(Integrator<ODE, odeintcpp::normalize>(      std::move(od), method, atol, rtol), inodes, tstates);
+      ll_res = calc_ll(Integrator<ODE, odeintcpp::normalize>(      std::move(od), method, atol, rtol), inodes, tstates, num_threads);
     } else {
-      ll_res = calc_ll(Integrator<ODE, odeintcpp::no_normalization>(std::move(od), method, atol, rtol), inodes, tstates);
+      ll_res = calc_ll(Integrator<ODE, odeintcpp::no_normalization>(std::move(od), method, atol, rtol), inodes, tstates, num_threads);
     }
 
 

src/secsse_loglik.h

@@ -101,28 +101,29 @@ namespace secsse {
   }
 
 
-  inline inodes_t<terse::inode_t> find_inte_nodes(const std::vector<phy_edge_t>& phy_edge, rvector<const int> ances, std::vector<std::vector<double>>& states) {
+  inline inodes_t<terse::inode_t> find_inte_nodes(const std::vector<phy_edge_t>& phy_edge, rvector<const int> ances, std::vector<std::vector<double>>& states, const int num_threads) {
     auto res = inodes_t<terse::inode_t>{ances.size()};
     auto comp = [](auto& edge, size_t val) { return edge.n < val; };
-    tbb::parallel_for<int>(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<int>(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 <typename RaIt>
   inline double normalize_loglik(RaIt first, RaIt last) {
     const auto sabs = std::accumulate(first, last, 0.0, [](const auto& s, const auto& x) {
@@ -230,16 +231,19 @@ namespace secsse {
   template <typename INTEGRATOR>
   inline calc_ll_res calc_ll(const INTEGRATOR& integrator,
                              inodes_t<terse::inode_t>& inodes,
-                             std::vector<std::vector<double>>& /* in/out */ states)
+                             std::vector<std::vector<double>>& /* in/out */ states,
+                             int num_threads)
   {
     const auto d = integrator.size();
     auto is_dirty = [](const auto& inode) {
       return inode.state->empty() && (inode.desc[0].state->empty() || inode.desc[1].state->empty());
     };
     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;
     }

src/tbb_stub.h

@@ -65,7 +65,10 @@ namespace tbb {
     }
   }
 
-
+  class task_arena {
+    task_arena(size_t /*num_threads*/) {});
+    template<typename Func> void execute(const Func& f) { f(); }
+  };
 } // namespce tbb
 
 

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 

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,23 +85,26 @@ 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,
                  traits = traits,
                  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,
                  traits = traits,
                  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)
 ```
 

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