Merge pull request #138 from rsetienne/develop

v4.0.5

.Rbuildignore

@@ -12,9 +12,6 @@
 ^pics$
 ^\.github$
 
-invoke_IW\.R
-test\.R
-test\.sh
 ^\.Rproj\.user$
 ^\.covrignore$
 ^LICENSE$

.github/workflows/R-CMD-check.yaml

@@ -2,17 +2,9 @@
 # https://community.rstudio.com/new-topic?category=Package%20development&tags=github-actions
 on:
   push:
-    branches:
-      - main
-      - master
-      - develop
-      - DAISIE_IW_B
-
   pull_request:
-    branches:
-      - main
-      - master
-      - develop
+  schedule:
+  - cron: "0 0 * * *"
 
 name: R-CMD-check
 
@@ -28,8 +20,8 @@ jobs:
         config:
           - {os: windows-latest, r: 'release'}
           - {os: macOS-latest, r: 'release'}
-          - {os: ubuntu-20.04, r: 'release', rspm: "https://packagemanager.rstudio.com/cran/__linux__/focal/latest"}
-          - {os: ubuntu-20.04, r: 'devel', rspm: "https://packagemanager.rstudio.com/cran/__linux__/focal/latest"}
+          - {os: ubuntu-latest, r: 'release', rspm: "https://packagemanager.rstudio.com/cran/__linux__/focal/latest"}
+          - {os: ubuntu-latest, r: 'devel', rspm: "https://packagemanager.rstudio.com/cran/__linux__/focal/latest"}
 
     env:
       R_REMOTES_NO_ERRORS_FROM_WARNINGS: true

.github/workflows/test-coverage.yaml

@@ -1,50 +1,27 @@
 on:
   push:
-    branches:
-      - main
-      - master
-      - develop
   pull_request:
-    branches:
-      - main
-      - master
-      - develop
+  schedule:
+  - cron: "0 0 * * *"
 
 name: test-coverage
 
 jobs:
   test-coverage:
-    runs-on: ubuntu-20.04
-
+    runs-on: ubuntu-latest
     env:
       GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }}
+
     steps:
       - uses: actions/checkout@v2
 
       - uses: r-lib/actions/setup-r@v1
-
-      - uses: r-lib/actions/setup-pandoc@v1
-
-      - name: Query dependencies
-        run: |
-          install.packages('remotes')
-          saveRDS(remotes::dev_package_deps(dependencies = TRUE), ".github/depends.Rds", version = 2)
-          writeLines(sprintf("R-%i.%i", getRversion()$major, getRversion()$minor), ".github/R-version")
-        shell: Rscript {0}
-
-      - name: Restore R package cache
-        uses: actions/cache@v2
         with:
-          path: ${{ env.R_LIBS_USER }}
-          key: ${{ runner.os }}-${{ hashFiles('.github/R-version') }}-1-${{ hashFiles('.github/depends.Rds') }}
-          restore-keys: ${{ runner.os }}-${{ hashFiles('.github/R-version') }}-1-
+          use-public-rspm: true
 
-      - name: Install dependencies
-        run: |
-          install.packages(c("remotes"))
-          remotes::install_deps(dependencies = TRUE)
-          remotes::install_cran("covr")
-        shell: Rscript {0}
+      - uses: r-lib/actions/setup-r-dependencies@v1
+        with:
+          extra-packages: covr
 
       - name: Test coverage
         run: covr::codecov()

DESCRIPTION

@@ -1,8 +1,8 @@
 Package: DAISIE
 Type: Package
 Title: Dynamical Assembly of Islands by Speciation, Immigration and Extinction
-Version: 4.0.2
-Date: 2021-05-26
+Version: 4.0.5
+Date: 2021-11-17
 Depends: R (>= 3.5.0)
 biocViews:
 SystemRequirements: C++14
@@ -115,4 +115,4 @@ Encoding: UTF-8
 VignetteBuilder: knitr
 URL: https://github.com/rsetienne/DAISIE
 BugReports: https://github.com/rsetienne/DAISIE/issues
-RoxygenNote: 7.1.1
+RoxygenNote: 7.1.2

NAMESPACE

@@ -2,6 +2,7 @@
 
 export(DAISIE_ExpEIN)
 export(DAISIE_IC)
+export(DAISIE_IW_num_threads)
 export(DAISIE_ML)
 export(DAISIE_ML_CS)
 export(DAISIE_ML_IW)
@@ -34,6 +35,7 @@ export(create_CS_version)
 export(create_area_pars)
 export(create_hyper_pars)
 export(create_trait_pars)
+export(daisie_odeint_cs)
 export(daisie_odeint_iw)
 import(Rcpp)
 importFrom(doParallel,registerDoParallel)

NEWS.md

@@ -1,16 +1,34 @@
+# DAISIE 4.0.5
+
+**N.B.: MacOS users may experience issues when installing DAISIE, especially when on MacOS Big Sur. If that is you case, please see [here](https://github.com/rsetienne/DAISIE/blob/6da0e3f65680d5f237345ef80935bda7541cf230/doc/DAISIE_macOS.md) for detailed installation instructions.**
+
+* CS model can be fit with C++ code, similarly to the IW model since version v3.2.0. To do so, set `methode` to:
+    * `"odeint::runge_kutta_cash_karp54"`
+    * `"odeint::runge_kutta_fehlberg78"` [default]
+    * `"odeint::runge_kutta_dopri5"`
+    * `"odeint::bulirsch_stoer"`
+
+without the `odeint::`-prefix, the `deSolve::ode` method is assumed. `deSolve::ode` using `"lsodes"` is still the default, as before. Make sure `CS_version = 0` (default) to use C++ ODEINT numerical integration.
+* Better handling of number of threads to use with `DAISIE_IW_num_threads()` for ML estimation of the IW model when more threads are available.
+* Fixes to numerical instability in conditional probability calculation in the IW model when the probability of extinction is very close to 1.
+* Introduce an approximation when the log conditional likelihood cannot be computed and add message when this is triggered.
+
 # DAISIE 4.0.2
+[![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.4813095.svg)](https://doi.org/10.5281/zenodo.4813095)
 
 **N.B.: MacOS users may experience issues when installing DAISIE, especially when on MacOS Big Sur. If that is you case, please see [here](https://github.com/rsetienne/DAISIE/blob/6da0e3f65680d5f237345ef80935bda7541cf230/doc/DAISIE_macOS.md) for detailed installation instructions.**
 
 * Suggest ggtree >= 3.0.0.
 
 # DAISIE 4.0.1
+[![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.4813094.svg)](https://doi.org/10.5281/zenodo.4813094)
 
 **N.B.: MacOS users may experience issues when installing DAISIE, especially when on MacOS Big Sur. If that is you case, please see [here](https://github.com/rsetienne/DAISIE/blob/6da0e3f65680d5f237345ef80935bda7541cf230/doc/DAISIE_macOS.md) for detailed installation instructions.**
 
 * Fix possibility of fitting CS model with IW likelihood on simulated data by setting `CS_version = 0`. Improve `CS_version` documentation.
 
 # DAISIE 4.0.0
+[![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.4813090.svg)](https://doi.org/10.5281/zenodo.4813090)
 
 **N.B.: MacOS users may experience issues when installing DAISIE, especially when on MacOS Big Sur. If that is you case, please see [here](https://github.com/rsetienne/DAISIE/blob/6da0e3f65680d5f237345ef80935bda7541cf230/doc/DAISIE_macOS.md) for detailed installation instructions.**
 

R/DAISIE-package.R

@@ -13,7 +13,9 @@
 #' \item Valente, L., Etienne, R. S., & Garcia-R, J. C. (2019). Deep macroevolutionary impact of humans on New Zealand's unique avifauna. Current Biology 29: 2563-2569. \doi{10.1016/j.cub.2019.06.058}.\cr
 #' \item Valente, L., Phillimore, A. B., Melo, M., Warren, B. H., Clegg, S. M., Havenstein, K., & Etienne, R. S. (2020). A simple dynamic model explains the diversity of island birds worldwide. Nature 579: 92-96. \doi{10.1038/s41586-020-2022-5}.\cr
 #' \item Hauffe, T., Delicado, D., Etienne, R.S., & Valente, L. (2020). Lake expansion elevates equilibrium diversity via increasing colonization. Journal of Biogeography 47: 1849–1860. \doi{10.1111/jbi.13914}.\cr
-#' } \cr
+#' \item Valente, L., Kristensen, N., Phillimore, A. B., & Etienne, R. S. (2021). Report of programming bugs in the DAISIE R package: consequences and correction. EcoEvoRxiv. \doi{10.32942/osf.io/w5ntf}.\cr
+#' \item Santos Neves, P., Lambert, J. W., Valente, L., & Etienne, R. S. (2021).The robustness of a simple dynamic model of island biodiversity to geological and eustatic change. bioRxiv 2021.07.26.453064. \doi{10.1101/2021.07.26.453064}.\cr
+#' }
 #' @keywords internal
 #' @import Rcpp
 #' @useDynLib DAISIE, .registration = TRUE

R/DAISIE_SR_loglik_CS.R

@@ -31,10 +31,8 @@ odeproc <- function(
                             method = method)
     } else
       if (times[1] < tshift & times[2] > tshift) {
-        #y = deSolve::ode(probs,c(times[1],tshift),fun,pars1,rtol = rtol,atol = atol,method = method)
-        y <- DAISIE_integrate(probs, c(times[1], tshift), fun, pars1, rtol = rtol, atol = atol, method = method)
-        probs <- y[2, 2:ncol(y)]
-        #y = deSolve::ode(probs,c(tshift,times[2]),fun,pars2,rtol = rtol,atol = atol,method = method)
+        probs <- DAISIE_integrate(probs, c(times[1], tshift), fun, pars1, rtol = rtol, atol = atol, method = method)
+
         y <- DAISIE_integrate(probs, c(tshift, times[2]), fun, pars2, rtol = rtol, atol = atol, method = method)
       }
   return(y)
@@ -147,8 +145,7 @@ DAISIE_SR_loglik_CS_M1 <- DAISIE_SR_loglik <- function(
       probs <- rep(0, 2 * lx + 1)
       probs[1] <- 1
       k1 <- 0
-      y <- odeproc(probs, brts[1:2], DAISIE_loglik_rhs, c(pars1, k1, ddep), rtol = reltolint, atol = abstolint, method = methode)
-      probs <- y[2, 2:(2 * lx + 2)]
+      probs <- odeproc(probs, brts[1:2], DAISIE_loglik_rhs, c(pars1, k1, ddep), rtol = reltolint, atol = abstolint, method = methode)
       cp <- checkprobs(lv = 2 * lx, loglik, probs, verbose); loglik <- cp[[1]]; probs <- cp[[2]]
       if (stac == 0) {
       # for stac = 0, the integration is from the origin of the island until the present
@@ -164,17 +161,15 @@ DAISIE_SR_loglik_CS_M1 <- DAISIE_SR_loglik <- function(
         # but there can be missing species
         # for stac = 5, we do exactly the same, but we evaluate the probability of an endemic species being present alone.
           probs[(lx + 1):(2 * lx)] <- 0
-          y <- odeproc(probs, brts[2:3], DAISIE_loglik_rhs, c(pars1, k1, ddep), rtol = reltolint, atol = abstolint, method = methode)
-          probs <- y[2, 2:(2 * lx + 2)]
+          probs <- odeproc(probs, brts[2:3], DAISIE_loglik_rhs, c(pars1, k1, ddep), rtol = reltolint, atol = abstolint, method = methode)
           cp <- checkprobs(lv = 2 * lx, loglik, probs, verbose); loglik <- cp[[1]]; probs <- cp[[2]]
           loglik <- loglik + log(probs[(stac == 1) * lx + (stac == 5) + 1 + missnumspec])
         } else {
         # for stac > 1, but not 5, integration is then from the colonization event until the first branching time (stac = 2 and 3) or the present (stac = 4). We add a set of equations for Q_M,n, the probability that the process is compatible with the data, and speciation has not happened; during this time immigration is not allowed because it would alter the colonization time. After speciation, colonization is allowed again (re-immigration)
         # all probabilities of states with the immigrant present are set to zero and all probabilities of states with endemics present are transported to the state with the colonist present waiting for speciation to happen. We also multiply by the (possibly diversity-dependent) immigration rate
           if (stac == 6 || stac == 7) {
             probs[(lx + 1):(2 * lx)] <- 0
-            y <- odeproc(probs, brts[2:3], DAISIE_loglik_rhs, c(pars1, k1, ddep), rtol = reltolint, atol = abstolint, method = methode)
-            probs <- y[2, 2:(2 * lx + 2)]
+            probs <- odeproc(probs, brts[2:3], DAISIE_loglik_rhs, c(pars1, k1, ddep), rtol = reltolint, atol = abstolint, method = methode)
             cp <- checkprobs(lv = 2 * lx, loglik, probs, verbose); loglik <- cp[[1]]; probs <- cp[[2]]
             k1 <- 1
           }
@@ -184,8 +179,7 @@ DAISIE_SR_loglik_CS_M1 <- DAISIE_SR_loglik <- function(
               gamvec[2:(lx + 1)] * probs[(lx + 1):(2 * lx)]
             probs[1:(2 * lx)] <- 0
             k1 <- 1
-            y <- odeproc(probs, brts[2:3], DAISIE_loglik_rhs2, c(pars1, k1, ddep), rtol = reltolint, atol = abstolint, method = methode)
-            probs <- y[2, 2:(3 * lx + 1)]
+            probs <- odeproc(probs, brts[2:3], DAISIE_loglik_rhs2, c(pars1, k1, ddep), rtol = reltolint, atol = abstolint, method = methode)
             cp <- checkprobs2(lx, loglik, probs, verbose); loglik <- cp[[1]]; probs <- cp[[2]]
           }
           if (stac == 4) {
@@ -215,8 +209,7 @@ DAISIE_SR_loglik_CS_M1 <- DAISIE_SR_loglik <- function(
               }
               for (k in startk:S1) {
                 k1 <- k - 1
-                y <- odeproc(probs, brts[k:(k + 1)], DAISIE_loglik_rhs, c(pars1, k1, ddep), rtol = reltolint, atol = abstolint, method = methode)
-                probs <- y[2, 2:(2 * lx + 2)]
+                probs <- odeproc(probs, brts[k:(k + 1)], DAISIE_loglik_rhs, c(pars1, k1, ddep), rtol = reltolint, atol = abstolint, method = methode)
                 cp <- checkprobs2(lx, loglik, probs, verbose); loglik <- cp[[1]]; probs <- cp[[2]]
                 if (k < S1) {
                   # speciation event