diff --git a/.github/workflows/R-CMD-check.yaml b/.github/workflows/R-CMD-check.yaml index fe099a77..a3d89ced 100644 --- a/.github/workflows/R-CMD-check.yaml +++ b/.github/workflows/R-CMD-check.yaml @@ -1,10 +1,8 @@ -# For help debugging build failures open an issue on the RStudio community with the 'github-actions' tag. -# https://community.rstudio.com/new-topic?category=Package%20development&tags=github-actions +# Workflow derived from https://github.com/r-lib/actions/tree/v2/examples +# Need help debugging build failures? Start at https://github.com/r-lib/actions#where-to-find-help on: push: pull_request: - schedule: - - cron: "0 0 * * *" name: R-CMD-check @@ -14,73 +12,37 @@ jobs: name: ${{ matrix.config.os }} (${{ matrix.config.r }}) + if: "contains(github.event.head_commit.message, '[run ci]') || (github.ref == 'refs/heads/master' || github.ref == 'refs/heads/develop')" + strategy: fail-fast: false matrix: config: + - {os: macos-latest, r: 'release'} - {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: 'devel', http-user-agent: 'release'} + - {os: ubuntu-latest, r: 'release'} env: - R_REMOTES_NO_ERRORS_FROM_WARNINGS: true - RSPM: ${{ matrix.config.rspm }} GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} + R_KEEP_PKG_SOURCE: yes steps: - - uses: actions/checkout@v2 + - uses: actions/checkout@v3 + + - uses: r-lib/actions/setup-pandoc@v2 - - uses: r-lib/actions/setup-r@v1 + - uses: r-lib/actions/setup-r@v2 with: r-version: ${{ matrix.config.r }} + http-user-agent: ${{ matrix.config.http-user-agent }} + use-public-rspm: true - - uses: r-lib/actions/setup-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: Cache R packages - if: runner.os != 'Windows' - uses: actions/cache@v2 + - uses: r-lib/actions/setup-r-dependencies@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- - - - name: Install system dependencies - if: runner.os == 'Linux' - run: | - while read -r cmd - do - eval sudo $cmd - done < <(Rscript -e 'writeLines(remotes::system_requirements("ubuntu", "20.04"))') - - - name: Install dependencies - run: | - remotes::install_deps(dependencies = TRUE) - remotes::install_cran("rcmdcheck") - if (.Platform$OS.type == "unix") remotes::install_cran("doMC") - shell: Rscript {0} - - - name: Check - env: - _R_CHECK_CRAN_INCOMING_REMOTE_: false - run: rcmdcheck::rcmdcheck(args = c("--no-manual", "--as-cran"), error_on = "warning", check_dir = "check") - shell: Rscript {0} - - - name: Show testthat output - if: always() - run: find check -name 'testthat.Rout*' -exec cat '{}' \; || true - shell: bash + extra-packages: any::rcmdcheck, doMC=?ignore + needs: check - - name: Upload check results - if: failure() - uses: actions/upload-artifact@main + - uses: r-lib/actions/check-r-package@v2 with: - name: ${{ runner.os }}-r${{ matrix.config.r }}-results - path: check + upload-snapshots: true diff --git a/.github/workflows/test-coverage.yaml b/.github/workflows/test-coverage.yaml index a869e0de..952a4f78 100644 --- a/.github/workflows/test-coverage.yaml +++ b/.github/workflows/test-coverage.yaml @@ -1,28 +1,50 @@ +# Workflow derived from https://github.com/r-lib/actions/tree/v2/examples +# Need help debugging build failures? Start at https://github.com/r-lib/actions#where-to-find-help on: push: pull_request: - schedule: - - cron: "0 0 * * *" name: test-coverage jobs: test-coverage: runs-on: ubuntu-latest + if: "contains(github.event.head_commit.message, '[run ci]') || (github.ref == 'refs/heads/master' || github.ref == 'refs/heads/develop')" + env: GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} steps: - - uses: actions/checkout@v2 + - uses: actions/checkout@v3 - - uses: r-lib/actions/setup-r@v1 + - uses: r-lib/actions/setup-r@v2 with: use-public-rspm: true - - uses: r-lib/actions/setup-r-dependencies@v1 + - uses: r-lib/actions/setup-r-dependencies@v2 with: - extra-packages: covr + extra-packages: any::covr + needs: coverage - name: Test coverage - run: covr::codecov() + run: | + covr::codecov( + quiet = FALSE, + clean = FALSE, + install_path = file.path(Sys.getenv("RUNNER_TEMP"), "package") + ) shell: Rscript {0} + + - name: Show testthat output + if: always() + run: | + ## -------------------------------------------------------------------- + find ${{ runner.temp }}/package -name 'testthat.Rout*' -exec cat '{}' \; || true + shell: bash + + - name: Upload test results + if: failure() + uses: actions/upload-artifact@v3 + with: + name: coverage-test-failures + path: ${{ runner.temp }}/package diff --git a/DESCRIPTION b/DESCRIPTION index b1773b26..28aab8fa 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,8 +1,8 @@ Package: DAISIE Type: Package Title: Dynamical Assembly of Islands by Speciation, Immigration and Extinction -Version: 4.2.1 -Date: 2022-06-08 +Version: 4.3.0 +Date: 2023-12-26 Depends: R (>= 3.5.0) biocViews: SystemRequirements: C++14 @@ -11,7 +11,7 @@ Imports: graphics, stats, utils, - DDD (>= 4.4), + DDD (>= 5.0), subplex, Matrix, tensor, @@ -21,11 +21,11 @@ Imports: doParallel, magrittr, parallel, - Rcpp (>= 1.0.5) + Rcpp (>= 1.0.10) LinkingTo: Rcpp, RcppEigen, - BH + BH (>= 1.81.0-1) Suggests: covr, testthat (>= 2.1.0), @@ -67,7 +67,7 @@ Authors@R: c( email = "j.w.l.lambert@rug.nl", comment = c(ORCID = "0000-0001-5218-3046")), person(given = "Pedro", - family = "Neves", + family = "Santos Neves", role = c("aut"), email = "p.m.santos.neves@rug.nl", comment = c(ORCID = "0000-0003-2561-4677")), @@ -84,7 +84,8 @@ Authors@R: c( person(given = "Hanno", family = "Hildenbrandt", email = "h.hildenbrandt@rug.nl", - role = c("aut")), + role = c("aut"), + comment = c(ORCID = "0000-0002-6784-1037")), person(given = "Torsten", family = "Hauffe", email = "torsten.hauffe@gmail.com", @@ -115,4 +116,4 @@ Encoding: UTF-8 VignetteBuilder: knitr URL: https://github.com/rsetienne/DAISIE BugReports: https://github.com/rsetienne/DAISIE/issues -RoxygenNote: 7.2.0 +RoxygenNote: 7.2.3 diff --git a/NAMESPACE b/NAMESPACE index 49285890..3cf5beee 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -14,6 +14,7 @@ export(DAISIE_SR_loglik_CS) export(DAISIE_SR_loglik_all) export(DAISIE_abm_factor) export(DAISIE_convertprobdist) +export(DAISIE_count_species) export(DAISIE_dataprep) export(DAISIE_loglik_CS) export(DAISIE_loglik_IW) @@ -32,12 +33,14 @@ export(DAISIE_sim_cr_shift) export(DAISIE_sim_relaxed_rate) export(DAISIE_sim_time_dep) export(DAISIE_sim_trait_dep) +export(DAISIE_sim_trait_dep_2K) export(are_area_pars) export(create_CS_version) export(create_area_pars) export(create_hyper_pars) export(create_pars) export(create_trait_pars) +export(create_trait_pars_2K) export(daisie_odeint_cs) export(daisie_odeint_iw) import(Rcpp) diff --git a/NEWS.md b/NEWS.md index 85b744f9..3edb86f7 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,16 @@ +# DAISIE 4.3.0 + +* Due to recent changes in CRAN policy that result in warnings due to the use of Rcpp related packages that link to libraries that use `sprintf` ('Rcpp' and 'BH'), 'DAISIE' now requires Rcpp v(>= 1.0.10) and BH v(>= 1.81.0-1). See https://github.com/RcppCore/Rcpp/pull/1236 and https://github.com/eddelbuettel/bh/pull/90 respectively. Also require package 'DDD' v(>= 5.0.0). +* Add code for the trait-dependent simulations considering two carrying capacities for each trait state. +* Improvements to the relaxed rate model fitting procedure. +* Add `DAISIE_count_species()` to count the number of species in a datalist or simulated data. +* Simply printing code by wrapping multiple instances into `print_parameters_and_loglik()`. +* Overall documentation improvements and updated references. +* Bug fix in `DAISIE_sim_cr_iw()`, which wrongly computed number of species at present. #147 +* Improve the likelihood calculation. Fixed incorrect likelihood present for type 2 case CS model inference if the values of lambda were high, in the C++ implementation. +* Overhaul GHA workflows to be more up to date and more conservative when to run. Feature branches require the tag +[run ci] be added to the commit message to run. + # DAISIE 4.2.1 * Bug fixes in `DAISIE_loglik_CS()`on the likelihood code for the continental sampling (probability of initial presence on the island). diff --git a/R/DAISIE-package.R b/R/DAISIE-package.R index 66572b99..95c1c772 100644 --- a/R/DAISIE-package.R +++ b/R/DAISIE-package.R @@ -14,9 +14,9 @@ #' \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 #' \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. \doi{10.1101/2021.07.26.453064}.\cr +#' \item Santos Neves, P., Lambert, J. W., Valente, L., & Etienne, R. S. (2022). The robustness of a simple dynamic model of island biodiversity to geological and sea-level change. Journal of Biogeography. \doi{10.1111/jbi.14519}.\cr #' \item Lambert, J. W., Santos Neves, P., Bilderbeek, R. L. C., Valente, L., Etienne, R. S. (2022). The effect of mainland dynamics on data and parameter estimates in island biogeography. bioRxiv. \doi{10.1101/2022.01.13.476210}.\cr -#' \item Xie, S., Valente, L., Etienne, R. S. (2022). A simple island biodiversity model is robust to trait dependence in diversification and colonization rates. biRrxiv. \doi{10.1101/2022.01.01.474685}.\cr +#' \item Xie, S., Valente, L., Etienne, R. S. (2023). Can we ignore trait-dependent colonization and diversification in island biogeography? Evolution. \doi{10.1093/evolut/qpad006}.\cr #' } #' @keywords internal #' @import Rcpp diff --git a/R/DAISIE_ML1.R b/R/DAISIE_ML1.R index 9eacd329..2f0df0d7 100644 --- a/R/DAISIE_ML1.R +++ b/R/DAISIE_ML1.R @@ -12,11 +12,13 @@ DAISIE_loglik_all_choosepar <- function(trparsopt, abstolint = 1E-16, reltolint = 1E-10) { all_no_shift <- 6:10 + non_oceanic_option <- FALSE if (max(idparsopt,-Inf) <= 6 && max(idparsfix,-Inf) <= 6 && (6 %in% idparsopt || 6 %in% idparsfix)) { idparsnoshift <- 7:11 all_no_shift <- 7:11 + non_oceanic_option <- TRUE } if (sum(idparsnoshift %in% (all_no_shift)) != 5) { trpars1 <- rep(0, 11) @@ -45,7 +47,7 @@ DAISIE_loglik_all_choosepar <- function(trparsopt, pars1[idparsnoshift] <- pars1[idparsnoshift - 5] } } - if (min(pars1) < 0) { + if (min(pars1) < 0 | (pars1[6] > 1 && non_oceanic_option == TRUE)) { loglik <- -Inf } else { loglik <- DAISIE_loglik_all( @@ -240,13 +242,19 @@ DAISIE_ML1 <- function( idpars <- sort(c(idparsopt, idparsfix, idparsnoshift, idparseq)) missnumspec <- unlist(lapply(datalist, function(list) {list$missing_species})) # nolint - if (sum(missnumspec) > (res - 1)) { + if (max(missnumspec) > (res - 1)) { cat( "The number of missing species is too large relative to the resolution of the ODE.\n") return(out2err) } + if (max(missnumspec) > res/10) { + warning( + "The number of missing species is quite low relative to the + resolution of the ODE.\n") + } + if ((length(idpars) != max(idpars))) { cat("The parameters to be optimized and/or fixed are incoherent.\n") return(out2err) @@ -391,20 +399,8 @@ DAISIE_ML1 <- function( df = length(initparsopt), conv = unlist(out$conv) ) - s1 <- sprintf( - "Maximum likelihood parameter estimates:\n lambda_c: %f\n mu: %f\n K: %f\n gamma: %f\n lambda_a: %f\n lambda_c2: %f\n mu2: %f\n K2: %f\n gamma2: %f\n lambda_a2: %f\n prop_type2: %f", - MLpars1[1], - MLpars1[2], - MLpars1[3], - MLpars1[4], - MLpars1[5], - MLpars1[6], - MLpars1[7], - MLpars1[8], - MLpars1[9], - MLpars1[10], - MLpars1[11] - ) + pars_to_print <- MLpars1[1:11] + parnames <- c('lambda^c','mu','K','gamma','lambda^a','lambda^c2','mu2','K2','gamma2','lambda^a2','prop_type2') } else if (all(all_no_shift == 7:11)) { out2 <- data.frame( lambda_c = MLpars1[1], @@ -417,15 +413,8 @@ DAISIE_ML1 <- function( df = length(initparsopt), conv = unlist(out$conv) ) - s1 <- sprintf( - "Maximum likelihood parameter estimates:\n lambda_c: %f\n mu: %f\n K: %f\n gamma: %f\n lambda_a: %f\n prob_init_pres: %f", - MLpars1[1], - MLpars1[2], - MLpars1[3], - MLpars1[4], - MLpars1[5], - MLpars1[6] - ) + pars_to_print <- MLpars1[1:6] + parnames <- c('lambda^c','mu','K','gamma','lambda^a','prob_init_pres') } else { out2 <- data.frame( lambda_c = MLpars1[1], @@ -437,17 +426,14 @@ DAISIE_ML1 <- function( df = length(initparsopt), conv = unlist(out$conv) ) - s1 <- sprintf( - "Maximum likelihood parameter estimates:\n lambda_c: %f\n mu: %f\n K: %f\n gamma: %f\n lambda_a: %f\n", - MLpars1[1], - MLpars1[2], - MLpars1[3], - MLpars1[4], - MLpars1[5] - ) + pars_to_print <- MLpars1[1:5] + parnames <- c('lambda^c','mu','K','gamma','lambda^a') } - s2 <- sprintf("Maximum loglikelihood: %f", ML) - cat("\n", s1, "\n", s2, "\n") + print_parameters_and_loglik(pars = pars_to_print, + loglik = ML, + verbose = TRUE, + parnames = parnames, + type = 'island_ML') if (eqmodel > 0) { M <- calcMN(datalist, MLpars1) ExpEIN <- DAISIE_ExpEIN(datalist[[1]]$island_age, MLpars1, M) # nolint start diff --git a/R/DAISIE_ML2.R b/R/DAISIE_ML2.R index 3cff782b..3bd3702c 100644 --- a/R/DAISIE_ML2.R +++ b/R/DAISIE_ML2.R @@ -196,9 +196,17 @@ DAISIE_ML2 <- function( MLpars1[i, 3] <- Inf } } - out2 <- data.frame(lambda_c = MLpars1[, 1], mu = MLpars1[, 2], K = MLpars1[, 3], gamma = MLpars1[, 4], lambda_a = MLpars1[, 5], loglik = ML, df = length(initparsopt), conv = unlist(out$conv)) - s1 <- sprintf("Maximum likelihood parameter estimates: %f", MLpars1) - s2 <- sprintf("Maximum loglikelihood: %f", ML) - cat("\n", s1, "\n", s2, "\n") + out2 <- data.frame(lambda_c = MLpars1[, 1], + mu = MLpars1[, 2], + K = MLpars1[, 3], + gamma = MLpars1[, 4], + lambda_a = MLpars1[, 5], + loglik = ML, + df = length(initparsopt), + conv = unlist(out$conv)) + print_parameters_and_loglik(pars = MLpars1, + loglik = ML, + verbose = TRUE, + type = 'multiple_island_ML') return(invisible(out2)) } diff --git a/R/DAISIE_ML4.R b/R/DAISIE_ML4.R index a08cbbe1..cc76e033 100644 --- a/R/DAISIE_ML4.R +++ b/R/DAISIE_ML4.R @@ -18,9 +18,14 @@ DAISIE_loglik_all_choosepar4 <- function(trparsopt, loglik <- -Inf } else { pars1 <- trpars1 / (1 - trpars1) - CS_version$sd <- pars1[6] + CS_version$par_sd <- pars1[6] pars1 <- pars1[-6] - if (min(pars1) < 0) { + pick <- which(c("cladogenesis", + "extinction", + "carrying_capacity", + "immigration", + "anagenesis") == CS_version$relaxed_par) + if (min(pars1) < 0 | pars1[pick] > CS_version$par_upper_bound) { loglik <- -Inf } else { loglik <- DAISIE::DAISIE_loglik_all( @@ -82,7 +87,9 @@ DAISIE_ML4 <- function( methode = "lsodes", optimmethod = "subplex", CS_version = create_CS_version(model = 2, - relaxed_par = "cladogenesis"), + relaxed_par = "cladogenesis", + par_sd = 0, + par_upper_bound = Inf), verbose = 0, tolint = c(1E-16, 1E-10), island_ontogeny = NA, @@ -229,17 +236,10 @@ DAISIE_ML4 <- function( df = length(initparsopt), conv = unlist(out$conv) ) - s1 <- sprintf( - "Maximum likelihood parameter estimates: lambda_c: %f, mu: %f, K: %f, - gamma: %f, lambda_a: %f, sd: %f", - MLpars1[1], - MLpars1[2], - MLpars1[3], - MLpars1[4], - MLpars1[5], - MLpars1[6] - ) - s2 <- sprintf("Maximum loglikelihood: %f", ML) - cat("\n", s1, "\n", s2, "\n") + print_parameters_and_loglik(pars = MLpars1[1:6], + loglik = ML, + verbose = TRUE, + parnames = c('lambda^c','mu','K','gamma','lambda^a','sd'), + type = 'island_ML') return(invisible(out2)) } diff --git a/R/DAISIE_ML_IW.R b/R/DAISIE_ML_IW.R index c130e12b..8fd289c4 100644 --- a/R/DAISIE_ML_IW.R +++ b/R/DAISIE_ML_IW.R @@ -171,8 +171,9 @@ DAISIE_ML_IW <- function( if (length(idparsfix) != 0) { MLpars1[idparsfix] <- parsfix } if (MLpars1[3] > 10 ^ 7){ MLpars1[3] <- Inf } out2 <- data.frame(lambda_c = MLpars1[1], mu = MLpars1[2], K = MLpars1[3], gamma = MLpars1[4], lambda_a = MLpars1[5], loglik = ML, df = length(initparsopt), conv = unlist(out$conv)) - s1 <- sprintf("Maximum likelihood parameter estimates: lambda_c: %f, mu: %f, K: %f, gamma: %f, lambda_a: %f", MLpars1[1], MLpars1[2], MLpars1[3], MLpars1[4], MLpars1[5]) - s2 <- sprintf("Maximum loglikelihood: %f", ML) - cat("\n", s1, "\n", s2, "\n") + print_parameters_and_loglik(pars = MLpars1[1:5], + loglik = ML, + verbose = TRUE, + type = 'island_ML') return(invisible(out2)) } diff --git a/R/DAISIE_MW_ML.R b/R/DAISIE_MW_ML.R index 3b2ae608..cf36cf04 100644 --- a/R/DAISIE_MW_ML.R +++ b/R/DAISIE_MW_ML.R @@ -477,89 +477,10 @@ DAISIE_MW_ML = function( MLpars1[idparsopt] = MLpars if(length(idparsfix) != 0) { MLpars1[idparsfix] = parsfix } if(MLpars1[5] > 10^7){ MLpars1[5] = Inf } - s1output <- function(MLpars1,distance_dep) - { - s1 <- switch(distance_dep, - power = sprintf('Maximum likelihood parameter estimates:\n - lambda_c = %f * A^ %f\n - mu = %f * A^ -%f\n - K = %f * A^ %f\n - M * gamma = %f * d^ -%f\n - lambda_a = %f * d^ %f\n',MLpars1[1],MLpars1[2],MLpars1[3],MLpars1[4],MLpars1[5],MLpars1[6],MLpars1[7],MLpars1[8],MLpars1[9],MLpars1[10]), - signoidal_col = sprintf('Maximum likelihood parameter estimates:\n - lambda_c = %f * A^ %f\n - mu = %f * A^ -%f\n - K = %f * A^ %f\n - M * gamma = %f * (1 - (d/%f)^%f / (1 + (d/%f)^%f )\n - lambda_a = %f * d^ %f\n',MLpars1[1],MLpars1[2],MLpars1[3],MLpars1[4],MLpars1[5],MLpars1[6],MLpars1[7],MLpars1[11],MLpars1[8],MLpars1[11],MLpars1[8],MLpars1[9],MLpars1[10]), - sigmoidal_ana = sprintf('Maximum likelihood parameter estimates:\n - lambda_c = %f * A^ %f\n - mu = %f * A^ -%f\n - K = %f * A^ %f\n - M * gamma = %f * d^ -%f\n - lambda_a = %f * (d/%f)^%f / (1 + (d/%f)^%f )\n',MLpars1[1],MLpars1[2],MLpars1[3],MLpars1[4],MLpars1[5],MLpars1[6],MLpars1[7],MLpars1[8],MLpars1[9],MLpars1[11],MLpars1[10],MLpars1[11],MLpars1[10]), - sigmoidal_clado = sprintf('Maximum likelihood parameter estimates:\n - lambda_c = %f * (d/%f)^%f / (1 + (d/%f)^%f )\n - mu = %f * A^ -%f\n - K = %f * A^ %f\n - M * gamma = %f * d^ -%f\n - lambda_a = %f * d^ %f\n',MLpars1[1],MLpars1[11],MLpars1[2],MLpars1[11],MLpars1[2],MLpars1[3],MLpars1[4],MLpars1[5],MLpars1[6],MLpars1[7],MLpars1[8],MLpars1[9],MLpars1[10]), - sigmoidal_col_ana = sprintf('Maximum likelihood parameter estimates:\n - lambda_c = %f * A^ %f\n - mu = %f * A^ -%f\n - K = %f * A^ %f\n - M * gamma = %f * (1 - (d/%f)^%f / (1 + (d/%f)^%f )\n - lambda_a = %f * (d/%f)^%f / (1 + (d/%f)^%f )\n',MLpars1[1],MLpars1[2],MLpars1[3],MLpars1[4],MLpars1[5],MLpars1[6],MLpars1[7],MLpars1[11],MLpars1[8],MLpars1[11],MLpars1[8],MLpars1[9],MLpars1[12],MLpars1[10],MLpars1[12],MLpars1[10]), - area_additive_clado = sprintf('Maximum likelihood parameter estimates:\n - lambda_c = %f * A^ %f * d^ %f \n - mu = %f * A^ -%f\n - K = %f * A^ %f\n - M * gamma = %f * d^ -%f\n - lambda_a = %f * d^ %f\n',MLpars1[1],MLpars1[2],MLpars1[11],MLpars1[3],MLpars1[4],MLpars1[5],MLpars1[6],MLpars1[7],MLpars1[8],MLpars1[9],MLpars1[10]), - area_interactive_clado = sprintf('Maximum likelihood parameter estimates:\n - lambda_c = %f * A^ (%f + %f * log(d)) \n - mu = %f * A^ -%f\n - K = %f * A^ %f\n - M * gamma = %f * d^ -%f\n - lambda_a = %f * d^ %f\n',MLpars1[1],MLpars1[2],MLpars1[11],MLpars1[3],MLpars1[4],MLpars1[5],MLpars1[6],MLpars1[7],MLpars1[8],MLpars1[9],MLpars1[10]), - area_interactive_clado0 = sprintf('Maximum likelihood parameter estimates:\n - lambda_c = %f * A^ (%f + %f * log(d)) \n - mu = %f * A^ -%f\n - K = %f * A^ %f\n - M * gamma = %f * d^ -%f\n - lambda_a = %f * d^ %f\n',MLpars1[1],MLpars1[2],MLpars1[11],MLpars1[3],MLpars1[4],MLpars1[5],MLpars1[6],MLpars1[7],MLpars1[8],MLpars1[9],MLpars1[10]), - area_interactive_clado1 = sprintf('Maximum likelihood parameter estimates:\n - lambda_c = %f * A^ (%f + d/%f)) \n - mu = %f * A^ -%f\n - K = %f * A^ %f\n - M * gamma = %f * d^ -%f\n - lambda_a = %f * d^ %f\n',MLpars1[1],MLpars1[2],MLpars1[11],MLpars1[3],MLpars1[4],MLpars1[5],MLpars1[6],MLpars1[7],MLpars1[8],MLpars1[9],MLpars1[10]), - area_interactive_clado2 = sprintf('Maximum likelihood parameter estimates:\n - lambda_c = %f * A^ (%f + d/(d + %f)) \n - mu = %f * A^ -%f\n - K = %f * A^ %f\n - M * gamma = %f * d^ -%f\n - lambda_a = %f * d^ %f\n',MLpars1[1],MLpars1[2],MLpars1[11],MLpars1[3],MLpars1[4],MLpars1[5],MLpars1[6],MLpars1[7],MLpars1[8],MLpars1[9],MLpars1[10]), - area_interactive_clado3 = sprintf('Maximum likelihood parameter estimates:\n - lambda_c = %f * (A + d/%f)^ %f\n \n - mu = %f * A^ -%f\n - K = %f * A^ %f\n - M * gamma = %f * d^ -%f\n - lambda_a = %f * d^ %f\n',MLpars1[1],MLpars1[11],MLpars1[2],MLpars1[3],MLpars1[4],MLpars1[5],MLpars1[6],MLpars1[7],MLpars1[8],MLpars1[9],MLpars1[10]), - area_interactive_clado4 = sprintf('Maximum likelihood parameter estimates:\n - lambda_c = %f * A^ (%f * d/(d + %f)) \n - mu = %f * A^ -%f\n - K = %f * A^ %f\n - M * gamma = %f * d^ -%f\n - lambda_a = %f * d^ %f\n',MLpars1[1],MLpars1[2],MLpars1[11],MLpars1[3],MLpars1[4],MLpars1[5],MLpars1[6],MLpars1[7],MLpars1[8],MLpars1[9],MLpars1[10]) - ) - return(s1) - } - s2 = sprintf('Maximum loglikelihood: %f',ML) if(is.element(distance_dep,distance_dep_options1)) { out2 = data.frame(lambda_c0 = MLpars1[1], y = MLpars1[2], mu_0 = MLpars1[3], x = MLpars1[4], K_0 = MLpars1[5], z = MLpars1[6], gamma_0 = MLpars1[7], alpha = MLpars1[8], lambda_a0 = MLpars1[9], beta = MLpars1[10], d_0 = MLpars1[11], loglik = ML, df = length(initparsopt), conv = unlist(out$conv)) - } else + } else { if(distance_dep == 'sigmoidal_col_ana') { out2 = data.frame(lambda_c0 = MLpars1[1], y = MLpars1[2], mu_0 = MLpars1[3], x = MLpars1[4], K_0 = MLpars1[5], z = MLpars1[6], gamma_0 = MLpars1[7], alpha = MLpars1[8], lambda_a0 = MLpars1[9], beta = MLpars1[10], d0_col = MLpars1[11], d0_ana = MLpars1[12], loglik = ML, df = length(initparsopt), conv = unlist(out$conv)) @@ -567,6 +488,11 @@ DAISIE_MW_ML = function( { out2 = data.frame(lambda_c0 = MLpars1[1], y = MLpars1[2], mu_0 = MLpars1[3], x = MLpars1[4], K_0 = MLpars1[5], z = MLpars1[6], gamma_0 = MLpars1[7], alpha = MLpars1[8], lambda_a0 = MLpars1[9], beta = MLpars1[10], loglik = ML, df = length(initparsopt), conv = unlist(out$conv)) } - cat("\n",s1output(MLpars1,distance_dep),"\n",s2,"\n") + } + print_parameters_and_loglik(pars = MLpars1, + loglik = ML, + verbose = TRUE, + type = 'global_ML', + distance_dep = distance_dep) return(invisible(out2)) } diff --git a/R/DAISIE_SR_loglik_CS.R b/R/DAISIE_SR_loglik_CS.R index adc596c0..e8a35ea7 100644 --- a/R/DAISIE_SR_loglik_CS.R +++ b/R/DAISIE_SR_loglik_CS.R @@ -231,15 +231,10 @@ DAISIE_SR_loglik_CS_M1 <- DAISIE_SR_loglik <- function( } } } - #print(head(probs,n = 15)) - if (pars2[4] >= 1) { - # if (length(pars1 > 5)) { - s1 <- sprintf("Status of colonist: %d, Parameters: %f %f %f %f %f ", stac, pars1[1], pars1[2], pars1[3], pars1[4], pars1[5]) - # } - s2 <- sprintf(", Loglikelihood: %f", loglik) - cat(s1, s2, "\n", sep = "") - utils::flush.console() - } + print_parameters_and_loglik(pars = c(stac,pars1[1:5]), + loglik = loglik, + verbose = pars2[4], + type = 'clade_loglik') return(as.numeric(loglik)) } diff --git a/R/DAISIE_count_species.R b/R/DAISIE_count_species.R new file mode 100644 index 00000000..b81c9e0a --- /dev/null +++ b/R/DAISIE_count_species.R @@ -0,0 +1,101 @@ +#' @title Count number of species in DAISIE datalist or simulated data. +#' +#' @description Calculates various island diversity metrics from island datasets. +#' +#' @inheritParams default_params_doc +#' @author Luis Valente +#' @seealso \code{\link{DAISIE_dataprep}}, +#' \code{\link{DAISIE_plot_island}} +#' +#' @return The output is a list containing the following items: +#' \item{clade_sizes_sorted}{ List showing the total number of species in each +#' island clade (including missing species). Each item [[i]] on the list +#' gives the sizes of all clades for a single island. If option +#' sort_clade_sizes = T, +#' the clade sizes for are +#' sorted by increasing number of species. If option sort_clade_sizes = F +#' the clade sizes are given in the same order as in the input datalist.} +#' \item{size_largest_clade}{ The total number of species in the largest +#' island clade +#' for each island.} +#' \item{mean_clade_size}{ Mean clade size (average of all island clades)} +#' \item{number_colonisations}{ The total number of colonisations (clades) on +#' each island.} +#' \item{total_number_species}{ The total number of species on each island. These +#' are the extant species at present, including missing species; in case of +#' simulations, this is the number of species present on the island at the +#' end of the +#' simulation.} +#' @examples +#' # Run function with clade sizes in the order they appear in the input data +#' data("NewZealand_birds_datalist") +#' species_count <- DAISIE_count_species(NewZealand_birds_datalist) +#' +#' # Run function with clade sizes in ascending order +#' species_count_sorted <- DAISIE_count_species( +#' NewZealand_birds_datalist, +#' sort_clade_sizes = TRUE +#' ) +#' @export +DAISIE_count_species <- function(islands, sort_clade_sizes = TRUE) { + if (length(grep("not_present", islands)) == 1) { + islands <- list(islands) + } + + replicates <- length(islands) + time <- islands[[1]][[1]]$island_age + + ###### Calculate overall species richness and + ## colonization statistics across all islands + number_colonists <- c() + number_species <- c() + size_largest_clade <- c() + mean_clade_size <- c() + clade_sizes <- c() + + for (i in 1:replicates) { + the_island <- islands[[i]] + ncols <- length(the_island) - 1 + number_colonists <- append(number_colonists, ncols) + + if (ncols == 0) { + number_species <- append(number_species, 0) + size_largest_clade <- append(size_largest_clade, 0) + mean_clade_size <- append(mean_clade_size, 0) + clade_sizes <- append(clade_sizes, 0) + } + + if (ncols > 0) { + btimes <- c() + miss_specs <- c() + + for (o in 2:length(the_island)) { + btimes <- append(btimes, length(the_island[[o]]$branching_times) - 1) + miss_specs <- append(miss_specs, the_island[[o]]$missing_species) + } + number_species <- append(number_species, sum(btimes, miss_specs)) + } + clade_sizes_dist <- btimes + miss_specs + size_largest_clade <- append(size_largest_clade, max(clade_sizes_dist)) + mean_clade_size <- append(mean_clade_size, round(mean(clade_sizes_dist), 2)) + if (isTRUE(sort_clade_sizes)) { + clade_sizes[i] <- list(sort(clade_sizes_dist)) + } else { + clade_sizes[i] <- list(clade_sizes_dist) + } + } + + overall_results <- list( + clade_sizes_sorted = clade_sizes, + size_largest_clade = size_largest_clade, + mean_clade_size = mean_clade_size, + number_colonisations = number_colonists, + total_number_species = number_species + ) + + if (sort_clade_sizes == F) { + names(overall_results)[1] <- "clade_sizes" + } + + return(overall_results) +} diff --git a/R/DAISIE_dataprep.R b/R/DAISIE_dataprep.R index fd5823d9..95a36060 100644 --- a/R/DAISIE_dataprep.R +++ b/R/DAISIE_dataprep.R @@ -186,13 +186,11 @@ DAISIE_dataprep = function(datatable, if(is.na(the_brts[1])){ the_brts<-island_age if(datatable[i,"Status"] == "Endemic" | datatable[i,"Status"] == "endemic" ){ - levels(datatable$Status) = append(levels(datatable$Status),"Endemic_MaxAge") datatable[i,"Status"] <-"Endemic_MaxAge"} if(datatable[i,"Status"] == "Non_endemic" | datatable[i,"Status"] == "Non_Endemic" | datatable[i,"Status"] == "NonEndemic" | datatable[i,"Status"] == "Nonendemic" | datatable[i,"Status"] == "nonendemic" | datatable[i,"Status"] == "non_endemic") { - levels(datatable$Status) = append(levels(datatable$Status),"Non_endemic_MaxAge") datatable[i,"Status"] <-"Non_endemic_MaxAge"} } @@ -205,13 +203,11 @@ DAISIE_dataprep = function(datatable, " is older than island age", sep = "")) } if(datatable[i,"Status"] == "Endemic" | datatable[i,"Status"] == "endemic" ){ - levels(datatable$Status) = append(levels(datatable$Status),"Endemic_MaxAge") datatable[i,"Status"] <-"Endemic_MaxAge"} if(datatable[i,"Status"] == "Non_endemic" | datatable[i,"Status"] == "Non_Endemic" | datatable[i,"Status"] == "NonEndemic" | datatable[i,"Status"] == "Nonendemic" | datatable[i,"Status"] == "nonendemic" | datatable[i,"Status"] == "non_endemic" ) { - levels(datatable$Status) = append(levels(datatable$Status),"Non_endemic_MaxAge") datatable[i,"Status"] <-"Non_endemic_MaxAge"} } diff --git a/R/DAISIE_loglik_CS.R b/R/DAISIE_loglik_CS.R index a108688a..d1c22a13 100644 --- a/R/DAISIE_loglik_CS.R +++ b/R/DAISIE_loglik_CS.R @@ -133,7 +133,6 @@ DAISIE_loglik_rhs <- function(t, x, parsvec) { -(muvec[il3 + 1] + lacvec[il3 + 1]) * nn[in3 + 1] * xx2[ix3] + -laavec[il3 + 1] * xx2[ix3] dx3 <- 0 - return(list(c(dx1,dx2,dx3))) } @@ -310,8 +309,9 @@ checkprobs2 <- function(lv, loglik, probs, verbose) { } else if (sum(probs) <= 0) { loglik <- -Inf } else { - loglik = loglik + log(sum(probs)) - probs = probs/sum(probs) + sp <- sum(sort(probs)) + loglik = loglik + log(sp) + probs = probs/sp } if (verbose) { message("Numerical issues encountered \n") @@ -495,11 +495,10 @@ DAISIE_loglik_CS_M1 <- DAISIE_loglik <- function(pars1, loglik <- -Inf return(loglik) } - N <- length(brts) - 1 - # exception for N = 1 in high lambda case lac <- pars1[1] - if(lac == Inf & missnumspec == 0 & length(pars1) == 5 & N > 1) { - loglik <- DAISIE_loglik_high_lambda(pars1, -brts, stac) + if(lac == Inf & missnumspec == 0 & length(pars1) == 5) { + if(verbose) warning('Infinite lambda detected') + loglik <- DAISIE_loglik_high_lambda(pars1, -brts, stac) } else { if (ddep == 1 | ddep == 11) { lx <- min( @@ -675,24 +674,16 @@ DAISIE_loglik_CS_M1 <- DAISIE_loglik <- function(pars1, } } - if(pars2[4] >= 1) - { - if (length(pars1) == 11) { # CHANGE - s1 <- sprintf('Status of colonist: %d, Parameters: %f %f %f %f %f %f', stac, pars1[5], pars1[6], pars1[7], pars1[8], pars1[9], pars1[10]) - } else { - s1 <- sprintf( - "Status of colonist: %d, Parameters: %f %f %f %f %f ", - stac, - pars1[1], - pars1[2], - pars1[3], - pars1[4], - pars1[5] - ) - } - s2 <- sprintf(", Loglikelihood: %f", loglik) - cat(s1, s2, "\n", sep = "") - utils::flush.console() + if (length(pars1) == 11) { # CHANGE + print_parameters_and_loglik(pars = c(stac, pars1[5:10]), # should this be 6:10, or 6:11? + loglik = loglik, + verbose = pars2[4], + type = 'clade_loglik') + } else { + print_parameters_and_loglik(pars = c(stac, pars1[1:5]), + loglik = loglik, + verbose = pars2[4], + type = 'clade_loglik') } if (is.na(loglik)) { message("NA in loglik encountered. Changing to -Inf.") @@ -936,7 +927,10 @@ DAISIE_loglik_CS <- DAISIE_loglik_all <- function( { message('Positive values of loglik encountered without possibility for approximation. Setting loglik to -Inf.') loglik <- -Inf - print_parameters_and_loglik(pars = pars, loglik = loglik, verbose = pars2[4]) + print_parameters_and_loglik(pars = pars, + loglik = loglik, + verbose = pars2[4], + type = 'island_loglik') return(loglik) } if (is.null(datalist[[1]]$not_present)) { @@ -1040,21 +1034,158 @@ DAISIE_loglik_CS <- DAISIE_loglik_all <- function( reltolint = reltolint) } } - print_parameters_and_loglik(pars = pars, loglik = loglik, verbose = pars2[4]) + print_parameters_and_loglik(pars = pars, + loglik = loglik, + verbose = pars2[4], + type = 'island_loglik') return(loglik) } -print_parameters_and_loglik <- function(pars, loglik, verbose) +print_parameters_and_loglik <- function(pars, + loglik, + verbose, + parnames = c("lambda^c", "mu", "K", "gamma", "lambda^a"), + type = 'island_loglik', + distance_dep = NULL) { if (verbose >= 1) { - s1 <- sprintf("Parameters: ") - s2 <- sprintf("%f ", pars) - s3 <- sprintf(", Loglikelihood: %f", loglik) - cat(s1, s2, s3, "\n", sep = "") + if(type == 'clade_loglik') { + s1a <- sprintf("Status of colonist: %d", pars[1]) + s1b <- sprintf("Parameters:") + s1 <- paste(s1a, s1b, sep = ', ') + s2 <- paste(sprintf("%f", pars[-1]), collapse = ', ') + s12 <- paste(s1, s2, collapse = ' ') + s3 <- paste(sprintf("Loglikelihood: %f", loglik), collapse = '') + cat(s12, s3, sep = ', ') + } else { + if(type == 'global_ML') { + s1 <- s1output(pars, distance_dep) + s3 <- sprintf("Maximum Loglikelihood: %f", loglik) + cat(s1,s3,sep = '\n') + } else { + if(is.null(ncol(pars))) { + lpars <- length(pars) + } else { + lpars <- ncol(pars) + } + if(lpars != length(parnames)) + { + warning('The vectors of parameters and parameter names have different lengths.') + parnames <- NULL + } else { + parnames <- paste(parnames, collapse = ', ') + } + if(type == 'island_ML') { + s1 <- sprintf("Maximum likelihood parameters: ") + s2 <- paste(sprintf("%f", pars), collapse = ', ') + s3 <- sprintf("Maximum Loglikelihood: %f", loglik) + cat(s1, parnames, s2, s3, sep = '\n') + } else { + if(type == 'multiple_island_ML') { + s1 <- sprintf("Maximum likelihood parameters: ") + s2 <- parnames + for(i in 1:nrow(pars)) { + s2 <- paste(s2,paste(sprintf("%f", pars[i,]), collapse = ', '), sep = '\n') + } + s3 <- sprintf("Maximum Loglikelihood: %f", loglik) + cat(s1, s2, s3, sep = '\n') + } else { + if(type == 'island_loglik') { + s1 <- sprintf("Parameters: ") + s2 <- paste(sprintf("%f", pars), collapse = ', ') + s3 <- sprintf("Loglikelihood: %f", loglik) + cat(s1, parnames, s2, s3, sep = '\n') + } else { + stop('Type of printing output unknown') + } + } + } + } + } + cat('\n') utils::flush.console() } } +s1output <- function(MLpars1,distance_dep) +{ + s1 <- switch(distance_dep, + power = sprintf('Maximum likelihood parameter estimates:\n + lambda_c = %f * A^ %f\n + mu = %f * A^ -%f\n + K = %f * A^ %f\n + M * gamma = %f * d^ -%f\n + lambda_a = %f * d^ %f\n',MLpars1[1],MLpars1[2],MLpars1[3],MLpars1[4],MLpars1[5],MLpars1[6],MLpars1[7],MLpars1[8],MLpars1[9],MLpars1[10]), + sigmoidal_col = sprintf('Maximum likelihood parameter estimates:\n + lambda_c = %f * A^ %f\n + mu = %f * A^ -%f\n + K = %f * A^ %f\n + M * gamma = %f * (1 - (d/%f)^%f / (1 + (d/%f)^%f )\n + lambda_a = %f * d^ %f\n',MLpars1[1],MLpars1[2],MLpars1[3],MLpars1[4],MLpars1[5],MLpars1[6],MLpars1[7],MLpars1[11],MLpars1[8],MLpars1[11],MLpars1[8],MLpars1[9],MLpars1[10]), + sigmoidal_ana = sprintf('Maximum likelihood parameter estimates:\n + lambda_c = %f * A^ %f\n + mu = %f * A^ -%f\n + K = %f * A^ %f\n + M * gamma = %f * d^ -%f\n + lambda_a = %f * (d/%f)^%f / (1 + (d/%f)^%f )\n',MLpars1[1],MLpars1[2],MLpars1[3],MLpars1[4],MLpars1[5],MLpars1[6],MLpars1[7],MLpars1[8],MLpars1[9],MLpars1[11],MLpars1[10],MLpars1[11],MLpars1[10]), + sigmoidal_clado = sprintf('Maximum likelihood parameter estimates:\n + lambda_c = %f * (d/%f)^%f / (1 + (d/%f)^%f )\n + mu = %f * A^ -%f\n + K = %f * A^ %f\n + M * gamma = %f * d^ -%f\n + lambda_a = %f * d^ %f\n',MLpars1[1],MLpars1[11],MLpars1[2],MLpars1[11],MLpars1[2],MLpars1[3],MLpars1[4],MLpars1[5],MLpars1[6],MLpars1[7],MLpars1[8],MLpars1[9],MLpars1[10]), + sigmoidal_col_ana = sprintf('Maximum likelihood parameter estimates:\n + lambda_c = %f * A^ %f\n + mu = %f * A^ -%f\n + K = %f * A^ %f\n + M * gamma = %f * (1 - (d/%f)^%f / (1 + (d/%f)^%f )\n + lambda_a = %f * (d/%f)^%f / (1 + (d/%f)^%f )\n',MLpars1[1],MLpars1[2],MLpars1[3],MLpars1[4],MLpars1[5],MLpars1[6],MLpars1[7],MLpars1[11],MLpars1[8],MLpars1[11],MLpars1[8],MLpars1[9],MLpars1[12],MLpars1[10],MLpars1[12],MLpars1[10]), + area_additive_clado = sprintf('Maximum likelihood parameter estimates:\n + lambda_c = %f * A^ %f * d^ %f \n + mu = %f * A^ -%f\n + K = %f * A^ %f\n + M * gamma = %f * d^ -%f\n + lambda_a = %f * d^ %f\n',MLpars1[1],MLpars1[2],MLpars1[11],MLpars1[3],MLpars1[4],MLpars1[5],MLpars1[6],MLpars1[7],MLpars1[8],MLpars1[9],MLpars1[10]), + area_interactive_clado = sprintf('Maximum likelihood parameter estimates:\n + lambda_c = %f * A^ (%f + %f * log(d)) \n + mu = %f * A^ -%f\n + K = %f * A^ %f\n + M * gamma = %f * d^ -%f\n + lambda_a = %f * d^ %f\n',MLpars1[1],MLpars1[2],MLpars1[11],MLpars1[3],MLpars1[4],MLpars1[5],MLpars1[6],MLpars1[7],MLpars1[8],MLpars1[9],MLpars1[10]), + area_interactive_clado0 = sprintf('Maximum likelihood parameter estimates:\n + lambda_c = %f * A^ (%f + %f * log(d)) \n + mu = %f * A^ -%f\n + K = %f * A^ %f\n + M * gamma = %f * d^ -%f\n + lambda_a = %f * d^ %f\n',MLpars1[1],MLpars1[2],MLpars1[11],MLpars1[3],MLpars1[4],MLpars1[5],MLpars1[6],MLpars1[7],MLpars1[8],MLpars1[9],MLpars1[10]), + area_interactive_clado1 = sprintf('Maximum likelihood parameter estimates:\n + lambda_c = %f * A^ (%f + d/%f)) \n + mu = %f * A^ -%f\n + K = %f * A^ %f\n + M * gamma = %f * d^ -%f\n + lambda_a = %f * d^ %f\n',MLpars1[1],MLpars1[2],MLpars1[11],MLpars1[3],MLpars1[4],MLpars1[5],MLpars1[6],MLpars1[7],MLpars1[8],MLpars1[9],MLpars1[10]), + area_interactive_clado2 = sprintf('Maximum likelihood parameter estimates:\n + lambda_c = %f * A^ (%f + d/(d + %f)) \n + mu = %f * A^ -%f\n + K = %f * A^ %f\n + M * gamma = %f * d^ -%f\n + lambda_a = %f * d^ %f\n',MLpars1[1],MLpars1[2],MLpars1[11],MLpars1[3],MLpars1[4],MLpars1[5],MLpars1[6],MLpars1[7],MLpars1[8],MLpars1[9],MLpars1[10]), + area_interactive_clado3 = sprintf('Maximum likelihood parameter estimates:\n + lambda_c = %f * (A + d/%f)^ %f\n \n + mu = %f * A^ -%f\n + K = %f * A^ %f\n + M * gamma = %f * d^ -%f\n + lambda_a = %f * d^ %f\n',MLpars1[1],MLpars1[11],MLpars1[2],MLpars1[3],MLpars1[4],MLpars1[5],MLpars1[6],MLpars1[7],MLpars1[8],MLpars1[9],MLpars1[10]), + area_interactive_clado4 = sprintf('Maximum likelihood parameter estimates:\n + lambda_c = %f * A^ (%f * d/(d + %f)) \n + mu = %f * A^ -%f\n + K = %f * A^ %f\n + M * gamma = %f * d^ -%f\n + lambda_a = %f * d^ %f\n',MLpars1[1],MLpars1[2],MLpars1[11],MLpars1[3],MLpars1[4],MLpars1[5],MLpars1[6],MLpars1[7],MLpars1[8],MLpars1[9],MLpars1[10]) + ) + return(s1) +} + DAISIE_integrate <- function(initprobs, tvec, rhs_func, @@ -1194,6 +1325,8 @@ DAISIE_ode_cs <- function( times = tvec, func = rhs_func, parms = parsvec, + atol = atol, + rtol = rtol, method = methode)[,1:(N + 1)] probs <- y[-1,-1] } else { diff --git a/R/DAISIE_loglik_IW.R b/R/DAISIE_loglik_IW.R index e55b5c94..72c7aa50 100644 --- a/R/DAISIE_loglik_IW.R +++ b/R/DAISIE_loglik_IW.R @@ -522,6 +522,9 @@ DAISIE_loglik_IW <- function( loglik <- loglik - logcond } } - print_parameters_and_loglik(pars = pars1, loglik = loglik, verbose = pars2[4]) + print_parameters_and_loglik(pars = pars1, + loglik = loglik, + verbose = pars2[4], + type = 'island_loglik') return(as.numeric(loglik)) } diff --git a/R/DAISIE_loglik_integrate.R b/R/DAISIE_loglik_integrate.R index 6e7cedaf..d8eb6baf 100644 --- a/R/DAISIE_loglik_integrate.R +++ b/R/DAISIE_loglik_integrate.R @@ -19,14 +19,14 @@ DAISIE_loglik_integrate <- function( verbose) { testit::assert(is.list(CS_version)) - par_sd <- CS_version$sd + par_sd <- CS_version$par_sd + par_upper_bound <- CS_version$par_upper_bound pick <- which(c("cladogenesis", "extinction", "carrying_capacity", "immigration", "anagenesis") == CS_version$relaxed_par) par_mean <- pars1[pick] - integrated_loglik <- integral_peak( logfun = Vectorize(DAISIE_loglik_integrand, vectorize.args = "DAISIE_par"), @@ -44,7 +44,12 @@ DAISIE_loglik_integrate <- function( verbose = verbose, pick = pick, par_mean = par_mean, - par_sd = par_sd) + par_sd = par_sd, + par_upper_bound = par_upper_bound) - + cum_rho(par_upper_bound = par_upper_bound, + DAISIE_dist_pars = list(par_mean = par_mean, + par_sd = par_sd) + ) return(integrated_loglik) } @@ -109,6 +114,27 @@ rho <- function(DAISIE_par, DAISIE_dist_pars) { return(gamma_den) } +#' Cumulative Gamma distribution parameterised with mean and standard deviation +#' +#' @inheritParams default_params_doc +#' +#' @return Numeric +#' @keywords internal +cum_rho <- function(par_upper_bound, DAISIE_dist_pars) { + + gamma_pars <- transform_gamma_pars( + par_mean = DAISIE_dist_pars$par_mean, + par_sd = DAISIE_dist_pars$par_sd) + + gamma_cum_prob <- stats::pgamma( + q = par_upper_bound, + shape = gamma_pars$shape, + scale = gamma_pars$scale, + log.p = TRUE) + + return(gamma_cum_prob) +} + #' @title Computes integral of a very peaked function, modified from the #' SADISA package #' @description computes the logarithm of the integral of exp(logfun) from 0 @@ -146,7 +172,8 @@ integral_peak <- function(logfun, verbose, pick, par_mean, - par_sd) { + par_sd, + par_upper_bound) { fun <- function(x) { exp(logfun(x, pars1, @@ -215,9 +242,10 @@ integral_peak <- function(logfun, par_sd = par_sd) if (gamma_pars$shape < 1) { lower <- min(exp(xmax), 1E-3) - pars1[pick] <- lower / 2 + pars1f <- pars1 + pars1f[pick] <- lower / 2 Q0 <- exp(DAISIE_loglik( - pars1 = pars1, + pars1 = pars1f, pars2 = pars2, brts = brts, stac = stac, @@ -237,17 +265,18 @@ integral_peak <- function(logfun, upper = exp(xmax), subdivisions = 1000, rel.tol = 1e-10, - abs.tol = 1e-10) - Q2 <- stats::integrate(f = fun, - lower = exp(xmax), - upper = Inf, - subdivisions = 1000, - rel.tol = 1e-10, - abs.tol = 1e-10) - Q1 <- Q1$value - Q2 <- Q2$value + abs.tol = 1e-10)$value + if(exp(xmax) < par_upper_bound) { + Q2 <- stats::integrate(f = fun, + lower = exp(xmax), + upper = par_upper_bound, + subdivisions = 1000, + rel.tol = 1e-10, + abs.tol = 1e-10)$value + } else { + Q2 <- 0 + } logQ <- log(Q0 + Q1 + Q2) - return(logQ) } diff --git a/R/DAISIE_rates.R b/R/DAISIE_rates.R index 4f44da80..77586314 100644 --- a/R/DAISIE_rates.R +++ b/R/DAISIE_rates.R @@ -156,7 +156,7 @@ update_rates_trait <- function(timeval, trait_pars = trait_pars, island_spec = island_spec ) - testit::assert(is.list(immig_rate)) + ext_rate <- get_ext_rate( mu = mu, hyper_pars = hyper_pars, @@ -166,14 +166,13 @@ update_rates_trait <- function(timeval, trait_pars = trait_pars, island_spec = island_spec ) - testit::assert(is.list(ext_rate)) + ana_rate <- get_ana_rate( laa = laa, num_immigrants = num_immigrants, trait_pars = trait_pars, island_spec = island_spec ) - testit::assert(is.list(ana_rate)) clado_rate <- get_clado_rate( lac = lac, hyper_pars = hyper_pars, @@ -183,19 +182,13 @@ update_rates_trait <- function(timeval, trait_pars = trait_pars, island_spec = island_spec ) - testit::assert(is.list(clado_rate)) - testit::assert(!is.null(trait_pars)) + + trans_rate <- get_trans_rate(trait_pars = trait_pars, island_spec = island_spec) - testit::assert(is.list(trans_rate)) - # trait_pars <- create_trait_pars(trans_rate = trans_rate$trans_rate1, - # immig_rate2 = immig_rate$immig_rate2, - # ext_rate2 = ext_rate$ext_rate2, - # ana_rate2 = ana_rate$ana_rate2, - # clado_rate2 = clado_rate$clado_rate2, - # trans_rate2 = trans_rate$trans_rate2, - # M2 = trait_pars$M2) + + rates <- list( immig_rate = immig_rate$immig_rate1, ext_rate = ext_rate$ext_rate1, @@ -390,6 +383,7 @@ get_ana_rate <- function(laa, intersect(which(island_spec[,4] == "I"), which(island_spec[,8] == "2")) ) + # testit::assert(is.numeric(ana_rate1)) # testit::assert(ana_rate1 >= 0) # testit::assert(is.numeric(ana_rate2)) @@ -462,28 +456,25 @@ get_clado_rate <- function(lac, # testit::assert(clado_rate >= 0) # testit::assert(is.numeric(clado_rate)) return(clado_rate) - } else { + }else{ num_spec_trait1 <- length(which(island_spec[, 8] == "1")) num_spec_trait2 <- length(which(island_spec[, 8] == "2")) - clado_rate1_pc <- max( - 0, lac * num_spec_trait1 * (1 - num_spec / K), - na.rm = TRUE) - clado_rate1_pc <- get_clado_rate_per_capita( - lac = lac, - d = d, - num_spec = num_spec, - K = K, - A = A - ) - clado_rate1 <- num_spec_trait1 * clado_rate1_pc - clado_rate2_pc <- get_clado_rate_per_capita( - lac = lac, - d = d, - num_spec = num_spec_trait2, - K = K, - A = A - ) - clado_rate2 <- num_spec_trait2* clado_rate2_pc + + if ("K2" %in% names(trait_pars)) { + clado_rate1 <- max( + 0, lac * num_spec_trait1 * (1 - num_spec_trait1 / K), + na.rm = TRUE) + clado_rate2 <- max( + 0, trait_pars$clado_rate2 * num_spec_trait2 * (1 - num_spec_trait2 / trait_pars$K2), + na.rm = TRUE) + } else { + clado_rate1 <- max( + 0, lac * num_spec_trait1 * (1 - num_spec / K), + na.rm = TRUE) + clado_rate2 <- max( + 0, trait_pars$clado_rate2 * num_spec_trait2 * (1 - num_spec / K), + na.rm = TRUE) + } # testit::assert(clado_rate1 >= 0) # testit::assert(clado_rate2 >= 0) # testit::assert(is.numeric(clado_rate1)) @@ -556,18 +547,20 @@ get_immig_rate <- function(gam, } else { mainland_n2 <- trait_pars$M2 gam2 <- trait_pars$immig_rate2 - immig_rate1 <- max(0, mainland_n * get_immig_rate_per_capita( - gam = gam, - num_spec = num_spec, - K = K, - A = A - ), na.rm = TRUE) - immig_rate2 <- max(0, mainland_n2 * get_immig_rate_per_capita( - gam = gam2, - num_spec = num_spec, - K = K, - A = A - ), na.rm = TRUE) + num_spec_trait1 <- length(which(island_spec[, 8] == "1")) + num_spec_trait2 <- length(which(island_spec[, 8] == "2")) + if ("K2" %in% names(trait_pars)) { + immig_rate1 <- max(c(mainland_n * gam * (1 - (num_spec_trait1 / K)), + 0), na.rm = TRUE) + immig_rate2 <- max(c(mainland_n2 * gam2 * (1 - (num_spec_trait2 / trait_pars$K2)), + 0), na.rm = TRUE) + } else { + immig_rate1 <- max(c(mainland_n * gam * (1 - (num_spec / K)), + 0), na.rm = TRUE) + immig_rate2 <- max(c(mainland_n2 * gam2 * (1 - (num_spec / K)), + 0), na.rm = TRUE) + } + # testit::assert(is.numeric(immig_rate1)) # testit::assert(immig_rate1 >= 0) # testit::assert(is.numeric(immig_rate2)) @@ -605,12 +598,20 @@ get_trans_rate <- function(trait_pars, num_spec_trait1 <- length(which(island_spec[, 8] == "1")) num_spec_trait2 <- length(which(island_spec[, 8] == "2")) } - trans_rate1 <- trait_pars$trans_rate * num_spec_trait1 - trans_rate2 <- trait_pars$trans_rate2 * num_spec_trait2 - testit::assert(is.numeric(trans_rate1)) - testit::assert(trans_rate1 >= 0) - testit::assert(is.numeric(trans_rate2)) - testit::assert(trans_rate2 >= 0) + if ("K2" %in% names(trait_pars)) { + trans_rate1 <- trait_pars$trans_rate * num_spec_trait1 * + (1 - (num_spec_trait2 / trait_pars$K2)) + trans_rate2 <- trait_pars$trans_rate2 * num_spec_trait2 * + (1 - (num_spec_trait1 / trait_pars$K)) + } else { + trans_rate1 <- trait_pars$trans_rate * num_spec_trait1 + trans_rate2 <- trait_pars$trans_rate2 * num_spec_trait2 + } + + # testit::assert(is.numeric(trans_rate1)) + # testit::assert(trans_rate1 >= 0) + # testit::assert(is.numeric(trans_rate2)) + # testit::assert(trans_rate2 >= 0) trans_list <- list(trans_rate1 = trans_rate1, trans_rate2 = trans_rate2) return(trans_list) diff --git a/R/DAISIE_sample_event_trait_dep.R b/R/DAISIE_sample_event_trait_dep.R index e3d381c0..e9c0c200 100644 --- a/R/DAISIE_sample_event_trait_dep.R +++ b/R/DAISIE_sample_event_trait_dep.R @@ -19,7 +19,7 @@ #' @author Shu Xie #' @keywords internal DAISIE_sample_event_trait_dep <- function(rates) { - testit::assert(are_rates(rates)) + # testit::assert(are_rates(rates)) possible_event <- sample(x = 1:10, size = 1, replace = FALSE, @@ -34,8 +34,8 @@ DAISIE_sample_event_trait_dep <- function(rates) { rates$clado_rate2, rates$trans_rate2) ) - testit::assert(is.numeric(possible_event)) - testit::assert(possible_event >= 1) + # testit::assert(is.numeric(possible_event)) + # testit::assert(possible_event >= 1) return(possible_event) } diff --git a/R/DAISIE_sim_core_trait_dep.R b/R/DAISIE_sim_core_trait_dep.R index a24b7094..d35e5c87 100644 --- a/R/DAISIE_sim_core_trait_dep.R +++ b/R/DAISIE_sim_core_trait_dep.R @@ -3,15 +3,15 @@ #' @inheritParams default_params_doc #' @keywords internal DAISIE_sim_core_trait_dep <- function( - time, - mainland_n, - pars, - island_ontogeny = 0, - sea_level = 0, - hyper_pars, - area_pars, - extcutoff = 1000, - trait_pars = NULL + time, + mainland_n, + pars, + island_ontogeny = 0, + sea_level = 0, + hyper_pars, + area_pars, + extcutoff = 1000, + trait_pars = NULL ) { #### Initialization #### @@ -31,16 +31,13 @@ DAISIE_sim_core_trait_dep <- function( colonisation is zero. Island cannot be colonised.") } - - #### what is the useage of maxspecID and how to set M1 and M2??#### - mainland_n2 <- trait_pars$M2 mainland_ntotal <- mainland_n + mainland_n2 testit::assert(mainland_ntotal > 0) if(mainland_n != 0){ mainland_spec <- seq(1, mainland_n, 1) }else{ - mainland_spec = c() + mainland_spec <- c() } maxspecID <- mainland_ntotal @@ -48,14 +45,6 @@ DAISIE_sim_core_trait_dep <- function( stt_table <- matrix(ncol = 7) colnames(stt_table) <- c("Time","nI","nA","nC","nI2","nA2","nC2") stt_table[1,] <- c(total_time,0,0,0,0,0,0) - # spec_tables <- list(stt_table = stt_table, - # init_nonend_spec = init_nonend_spec, - # init_end_spec = init_end_spec, - # mainland_spec = mainland_spec, - # island_spec = island_spec) - # stt_table <- spec_tables$stt_table - # mainland_spec <- spec_tables$mainland_spec - # island_spec <- spec_tables$island_spec lac <- pars[1] mu <- pars[2] K <- pars[3] @@ -63,13 +52,7 @@ DAISIE_sim_core_trait_dep <- function( laa <- pars[5] num_spec <- length(island_spec[, 1]) - num_spec_trait1 <- length(which(island_spec[,8] == "1")) - num_spec_trait2 <- length(which(island_spec[,8] == "2")) num_immigrants <- length(which(island_spec[, 4] == "I")) - num_immigrants_trait1 <- length(intersect(which(island_spec[, 4] == "I"), - which(island_spec[, 8] == "1"))) - num_immigrants_trait2 <- length(intersect(which(island_spec[, 4] == "I"), - which(island_spec[, 8] == "2"))) #### Start Monte Carlo iterations #### while (timeval < total_time) { @@ -92,7 +75,6 @@ DAISIE_sim_core_trait_dep <- function( island_spec = island_spec, trait_pars = trait_pars ) - testit::assert(are_rates(rates)) timeval_and_dt <- calc_next_timeval( max_rates = rates, timeval = timeval @@ -100,26 +82,6 @@ DAISIE_sim_core_trait_dep <- function( timeval <- timeval_and_dt$timeval if (timeval < total_time) { - rates <- update_rates( - timeval = timeval, - total_time = total_time, - gam = gam, - laa = laa, - lac = lac, - mu = mu, - hyper_pars = hyper_pars, - area_pars = area_pars, - K = K, - num_spec = num_spec, - num_immigrants = num_immigrants, - mainland_n = mainland_n, - extcutoff = NULL, - island_ontogeny = 0, - sea_level = 0, - island_spec = island_spec, - trait_pars = trait_pars - ) - testit::assert(are_rates(rates)) possible_event <- DAISIE_sample_event_trait_dep( rates = rates ) @@ -161,7 +123,7 @@ DAISIE_sim_core_trait_dep <- function( island_spec = island_spec, mainland_n = mainland_n, trait_pars = trait_pars) - ordered_stt_times <- sort(island$stt_table[, 1], decreasing = TRUE) - testit::assert(all(ordered_stt_times == island$stt_table[, 1])) + # ordered_stt_times <- sort(island$stt_table[, 1], decreasing = TRUE) + # testit::assert(all(ordered_stt_times == island$stt_table[, 1])) return(island) } diff --git a/R/DAISIE_sim_cr_iw.R b/R/DAISIE_sim_cr_iw.R index 8090f2f8..d2db9bad 100644 --- a/R/DAISIE_sim_cr_iw.R +++ b/R/DAISIE_sim_cr_iw.R @@ -30,7 +30,8 @@ DAISIE_sim_cr_iw <- function(total_time, hyper_pars = hyper_pars, area_pars = area_pars ) - stac_vec <- unlist(island_replicates)[which(names(unlist(island_replicates)) == "taxon_list.stac")] + temp_island_replicates <- island_replicates[[rep]] + stac_vec <- unlist(temp_island_replicates)[which(names(unlist(temp_island_replicates)) == "taxon_list.stac")] present <- which(stac_vec != 0) number_present <- length(present) } diff --git a/R/DAISIE_sim_relaxed_rate.R b/R/DAISIE_sim_relaxed_rate.R index 92cc9ddd..1c2271df 100644 --- a/R/DAISIE_sim_relaxed_rate.R +++ b/R/DAISIE_sim_relaxed_rate.R @@ -76,8 +76,8 @@ #' carr_cap <- Inf #' immig_rate <- 0.005 #' ana_rate <- 1 -#' sd <- 1 -#' sim_pars <- c(clado_rate, ext_rate, carr_cap, immig_rate, ana_rate, sd) +#' par_sd <- 1 +#' sim_pars <- c(clado_rate, ext_rate, carr_cap, immig_rate, ana_rate, par_sd) #' set.seed(1) #' island_replicates <- DAISIE_sim_relaxed_rate( #' time = 1, diff --git a/R/DAISIE_sim_trait_dep.R b/R/DAISIE_sim_trait_dep.R index 3a217615..4beeb503 100644 --- a/R/DAISIE_sim_trait_dep.R +++ b/R/DAISIE_sim_trait_dep.R @@ -1,53 +1,41 @@ -#' @title Simulate islands with given parameters. +#' @title Simulate islands with given trait-dependent parameters. #' @description This function simulates islands with given cladogenesis, -#' extinction, Kprime, immigration and anagenesis parameters. If a single -#' parameter set is provided (5 parameters) it simulates islands where all -#' species have the same macro-evolutionary process. If two paramater sets -#' (10 parameters) are provided, it simulates islands where two different -#' macro-evolutionary processes operate, one applying to type 1 species and -#' other to type 2 species. If two parameter sets and a time shift (11 -#' parameters) are provided, it simulates islands where at the given time -#' a shift between the parameter sets will occur. +#' extinction, K, immigration and anagenesis parameters for binary states. #' #' Returns R list object that contains the simulated islands #' #' @inheritParams default_params_doc #' -#' @return Each simulated dataset is an element of the list, which can be -#' called using [[x]]. For example if the object is called island_replicates, -#' the last replicates is a list in itself. The first (e.g. -#' \code{island_replicates[[x]][[1]]}) element of that list has the following -#' components: \cr \code{$island_age} - the island age \cr Then, depending on -#' whether a distinction between types is made, we have:\cr \code{$not_present} -#' - the number of mainland lineages that are not present on the island \cr -#' or:\cr \code{$not_present_type1} - the number of mainland lineages of type 1 -#' that are not present on the island \cr \code{$not_present_type2} - the -#' number of mainland lineages of type 2 that are not present on the island \cr -#' \code{$stt_all} - STT table for all species on the island (nI - number of -#' non-endemic species; nA - number of anagenetic species, nC - number of -#' cladogenetic species, present - number of independent colonisations present -#' )\cr \code{$stt_stt_type1} - STT table for type 1 species on the island - -#' only if 2 types of species were simulated (nI - number of non-endemic -#' species; nA - number of anagenetic species, nC - number of cladogenetic -#' species, present - number of independent colonisations present )\cr -#' \code{$stt_stt_type2} - STT table for type 2 species on the island - only if -#' 2 types of species were simulated (nI - number of non-endemic species; nA - -#' number of anagenetic species, nC - number of cladogenetic species, present - -#' number of independent colonisations present )\cr \code{$brts_table} - Only -#' for simulations under 'IW'. Table containing information on order of events -#' in the data, for use in maximum likelihood optimization.)\cr -#' -#' The subsequent elements of the list each contain information on a single -#' colonist lineage on the island and has 4 components:\cr -#' \code{$branching_times} - island age and stem age of the population/species -#' in the case of Non-endemic, Non-endemic_MaxAge and Endemic anagenetic -#' species. For cladogenetic species these should be island age and branching -#' times of the radiation including the stem age of the radiation.\cr -#' \code{$stac} - the status of the colonist \cr * Non_endemic_MaxAge: 1 \cr * -#' ndemic: 2 \cr * Endemic&Non_Endemic: 3 \cr * Non_endemic: 4 \cr -#' \code{$missing_species} - number of island species that were not sampled for -#' particular clade (only applicable for endemic clades) \cr \code{$type_1or2} -#' - whether the colonist belongs to type 1 or type 2 \cr +#' @return +#' A list. The highest level of the least corresponds to each individual +#' replciate. The first element of each replicate is composed of island +#' information containing: +#' \itemize{ +#' \item{\code{$island_age}: A numeric with the island age.} +#' \item{\code{$not_present}: A numeric with the number of mainland lineages +#' that are not present on the island.} +#' \item{\code{$stt_all}: STT table for all species on the island +#' (nI - number of non-endemic species; nA - number of anagenetic species, +#' nC - number of cladogenetic species, present - number of independent +#' colonisations present)} +#' \item{\code{$brts_table}: Only for simulations under \code{"IW"}. Table +#' containing information on order of events in the data, for use in maximum +#' likelihood optimization.).} +#' } +#' The subsequent elements of the list pertaining to each replcate contain +#' information on a single colonist lineage on the island and have 4 components: +#' \itemize{ +#' \item{\code{$branching_times}: island age and stem age of the +#' population/species in the case of Non-endemic, Non-endemic_MaxAge and +#' Endemic anagenetic species. For cladogenetic species these should +#' be island age and branching times of the radiation including the +#' stem age of the radiation.} +#' \item{\code{$stac}: An integer ranging from 1 to 4 +#' indicating the status of the colonist:} +#' \item{\code{$missing_species}: number of island species that were +#' not sampled for particular clade (only applicable for endemic clades)} +#' \item{\code{$type_1or2}: whether the colonist belongs to type 1 or type 2} +#' } #' @author Luis Valente and Albert Phillimore #' @seealso \code{\link{DAISIE_format_CS}} \code{\link{DAISIE_plot_sims}} #' @references Valente, L.M., A.B. Phillimore and R.S. Etienne (2015). @@ -64,7 +52,6 @@ DAISIE_sim_trait_dep <- function( pars, replicates, divdepmodel = "CS", - num_guilds = NULL, sample_freq = 25, plot_sims = TRUE, island_ontogeny = "const", @@ -84,20 +71,6 @@ DAISIE_sim_trait_dep <- function( trait_pars = NULL, ... ) { - testit::assert( - "island_ontogeny is not valid input. Specify 'const',\n - or 'beta'", is_island_ontogeny_input(island_ontogeny) - ) - testit::assert( - "sea_level is not valid input. Specify 'const, \n or 'sine'", - is_sea_level_input(sea_level) - ) - - testit::assert( - "length(pars) is not five", - length(pars) == 5 - ) - total_time <- time island_replicates <- list() island_ontogeny <- translate_island_ontogeny(island_ontogeny) @@ -141,10 +114,11 @@ DAISIE_sim_trait_dep <- function( number_present <- 0 } while (number_present < cond) { - if(M == 0){ - if(is.null(trait_pars)){ - stop("There is no species on mainland.") - }else{ ## only have state2 species on mainland + if(M == 0 || is.null(trait_pars)){ + stop("One state exist on mainland, should use constant rate DAISIE.") + }else{ + for (m_spec in 1:M) { + ### M1 = 1, M2 = 0 trait_pars_onecolonize <- create_trait_pars( trans_rate = trait_pars$trans_rate, immig_rate2 = trait_pars$immig_rate2, @@ -152,31 +126,7 @@ DAISIE_sim_trait_dep <- function( ana_rate2 = trait_pars$ana_rate2, clado_rate2 = trait_pars$clado_rate2, trans_rate2 = trait_pars$trans_rate2, - M2 = 1) - for (m_spec in 1:trait_pars$M2) { - full_list[[m_spec]] <- DAISIE_sim_core_trait_dep( - time = total_time, - mainland_n = 0, - pars = pars, - island_ontogeny = island_ontogeny, - sea_level = sea_level, - hyper_pars = hyper_pars, - area_pars = area_pars, - extcutoff = extcutoff, - trait_pars = trait_pars_onecolonize - ) - } - } - }else{ - trait_pars_addcol <- create_trait_pars( - trans_rate = trait_pars$trans_rate, - immig_rate2 = trait_pars$immig_rate2, - ext_rate2 = trait_pars$ext_rate2, - ana_rate2 = trait_pars$ana_rate2, - clado_rate2 = trait_pars$clado_rate2, - trans_rate2 = trait_pars$trans_rate2, - M2 = 0) - for (m_spec in 1:M) { + M2 = 0) full_list[[m_spec]] <- DAISIE_sim_core_trait_dep( time = total_time, mainland_n = 1, @@ -186,11 +136,11 @@ DAISIE_sim_trait_dep <- function( hyper_pars = hyper_pars, area_pars = area_pars, extcutoff = extcutoff, - trait_pars = trait_pars_addcol + trait_pars = trait_pars_onecolonize ) } - for(m_spec in (M + 1):(M + trait_pars$M2)) - { + for(m_spec in (M + 1):(M + trait_pars$M2)) { + ### M1 = 0, M2 = 1 trait_pars_onecolonize <- create_trait_pars( trans_rate = trait_pars$trans_rate, immig_rate2 = trait_pars$immig_rate2, @@ -231,42 +181,6 @@ DAISIE_sim_trait_dep <- function( ) } - #### GW #### - if (divdepmodel == "GW") { - if (!is.numeric(num_guilds)) { - stop("num_guilds must be numeric") - } - guild_size <- M / num_guilds - testit::assert(num_guilds < M) - testit::assert(M %% num_guilds == 0) - for (rep in 1:replicates) { - island_replicates[[rep]] <- list() - full_list <- list() - for (m_spec in 1:num_guilds) { - full_list[[m_spec]] <- DAISIE_sim_core_trait_dep( - time = total_time, - mainland_n = guild_size, - pars = pars, - island_ontogeny = island_ontogeny, - sea_level = sea_level, - hyper_pars = hyper_pars, - area_pars = area_pars, - extcutoff = extcutoff - ) - } - island_replicates[[rep]] <- full_list - if (verbose == TRUE) { - print(paste("Island replicate ", rep, sep = "")) - } - } - island_replicates <- DAISIE_format_GW(island_replicates = island_replicates, - time = total_time, - M = M, - sample_freq = sample_freq, - num_guilds = num_guilds, - verbose = verbose) - } - if (plot_sims == TRUE) { DAISIE_plot_sims( island_replicates = island_replicates, diff --git a/R/DAISIE_sim_trait_dep_2K.R b/R/DAISIE_sim_trait_dep_2K.R new file mode 100644 index 00000000..5496884e --- /dev/null +++ b/R/DAISIE_sim_trait_dep_2K.R @@ -0,0 +1,195 @@ +#' @title Simulate islands with given trait-dependent parameters. +#' @description This function simulates islands with given cladogenesis, +#' extinction, K, immigration and anagenesis parameters. In this version, +#' rates and K are both trait-dependent. +#' +#' Returns R list object that contains the simulated islands +#' +#' @inheritParams default_params_doc +#' +#' @return +#' A list. The highest level of the least corresponds to each individual +#' replciate. The first element of each replicate is composed of island +#' information containing: +#' \itemize{ +#' \item{\code{$island_age}: A numeric with the island age.} +#' \item{\code{$not_present}: A numeric with the number of mainland lineages +#' that are not present on the island.} +#' \item{\code{$stt_all}: STT table for all species on the island +#' (nI - number of non-endemic species; nA - number of anagenetic species, +#' nC - number of cladogenetic species, present - number of independent +#' colonisations present)} +#' \item{\code{$brts_table}: Only for simulations under \code{"IW"}. Table +#' containing information on order of events in the data, for use in maximum +#' likelihood optimization.).} +#' } +#' The subsequent elements of the list pertaining to each replcate contain +#' information on a single colonist lineage on the island and have 4 components: +#' \itemize{ +#' \item{\code{$branching_times}: island age and stem age of the +#' population/species in the case of Non-endemic, Non-endemic_MaxAge and +#' Endemic anagenetic species. For cladogenetic species these should +#' be island age and branching times of the radiation including the +#' stem age of the radiation.} +#' \item{\code{$stac}: An integer ranging from 1 to 4 +#' indicating the status of the colonist:} +#' \item{\code{$missing_species}: number of island species that were +#' not sampled for particular clade (only applicable for endemic clades)} +#' \item{\code{$type_1or2}: whether the colonist belongs to type 1 or type 2} +#' } +#' @author Luis Valente and Albert Phillimore +#' @seealso \code{\link{DAISIE_format_CS}} \code{\link{DAISIE_plot_sims}} +#' @references Valente, L.M., A.B. Phillimore and R.S. Etienne (2015). +#' Equilibrium and non-equilibrium dynamics simultaneously operate in the +#' Galapagos islands. Ecology Letters 18: 844-852. +#' Hauffe, T., D. Delicado, R.S. Etienne and L. Valente (submitted). +#' Lake expansion increases equilibrium diversity via the target effect of +#' island biogeography. +#' @keywords models +#' @export +DAISIE_sim_trait_dep_2K <- function( + time, + M, + pars, + replicates, + divdepmodel = "CS", + sample_freq = 25, + plot_sims = TRUE, + island_ontogeny = "const", + sea_level = "const", + hyper_pars = create_hyper_pars(d = 0, x = 0), + area_pars = DAISIE::create_area_pars( + max_area = 1, + current_area = 1, + proportional_peak_t = 0, + total_island_age = 0, + sea_level_amplitude = 0, + sea_level_frequency = 0, + island_gradient_angle = 0), + extcutoff = 1000, + cond = 0, + verbose = TRUE, + trait_pars = NULL, + ... +) { + total_time <- time + island_replicates <- list() + island_ontogeny <- translate_island_ontogeny(island_ontogeny) + sea_level <- translate_sea_level(sea_level) + + #### IW #### + if (divdepmodel == "IW") { + for (rep in 1:replicates) { + island_replicates[[rep]] <- DAISIE_sim_core_trait_dep( + time = total_time, + mainland_n = M, + pars = pars, + island_ontogeny = island_ontogeny, + sea_level = sea_level, + hyper_pars = hyper_pars, + area_pars = area_pars, + extcutoff = extcutoff, + trait_pars = trait_pars + ) + if (verbose == TRUE) { + print(paste("Island replicate ", rep, sep = "")) + } + } + island_replicates <- DAISIE_format_IW( + island_replicates = island_replicates, + time = total_time, + M = M, + sample_freq = sample_freq, + verbose = verbose, + trait_pars = trait_pars) + } + + #### CS #### + if (divdepmodel == "CS") { + for (rep in 1:replicates) { + island_replicates[[rep]] <- list() + full_list <- list() + if (cond == 0) { + number_present <- -1 + } else { + number_present <- 0 + } + while (number_present < cond) { + if(M == 0 || is.null(trait_pars)){ + stop("One state exist on mainland, should use constant rate DAISIE.") + }else{ + for (m_spec in 1:M) { + ### M1 = 1, M2 = 0 + trait_pars_onecolonize <- create_trait_pars_2K( + trans_rate = trait_pars$trans_rate, + immig_rate2 = trait_pars$immig_rate2, + ext_rate2 = trait_pars$ext_rate2, + ana_rate2 = trait_pars$ana_rate2, + clado_rate2 = trait_pars$clado_rate2, + trans_rate2 = trait_pars$trans_rate2, + M2 = 0, + K2 = trait_pars$K2) + full_list[[m_spec]] <- DAISIE_sim_core_trait_dep( + time = total_time, + mainland_n = 1, + pars = pars, + island_ontogeny = island_ontogeny, + sea_level = sea_level, + hyper_pars = hyper_pars, + area_pars = area_pars, + extcutoff = extcutoff, + trait_pars = trait_pars_onecolonize + ) + } + for(m_spec in (M + 1):(M + trait_pars$M2)) { + ### M1 = 0, M2 = 1 + trait_pars_onecolonize <- create_trait_pars_2K( + trans_rate = trait_pars$trans_rate, + immig_rate2 = trait_pars$immig_rate2, + ext_rate2 = trait_pars$ext_rate2, + ana_rate2 = trait_pars$ana_rate2, + clado_rate2 = trait_pars$clado_rate2, + trans_rate2 = trait_pars$trans_rate2, + M2 = 1, + K2 = trait_pars$K2) + full_list[[m_spec]] <- DAISIE_sim_core_trait_dep( + time = total_time, + mainland_n = 0, + pars = pars, + island_ontogeny = island_ontogeny, + sea_level = sea_level, + hyper_pars = hyper_pars, + area_pars = area_pars, + extcutoff = extcutoff, + trait_pars = trait_pars_onecolonize + ) + } + } + stac_vec <- unlist(full_list)[which(names(unlist(full_list)) == "stac")] + present <- which(stac_vec != 0) + number_present <- length(present) + } + island_replicates[[rep]] <- full_list + if (verbose == TRUE) { + print(paste("Island replicate ", rep, sep = "")) + } + } + island_replicates <- DAISIE_format_CS( + island_replicates = island_replicates, + time = total_time, + M = M, + sample_freq = sample_freq, + verbose = verbose, + trait_pars = trait_pars + ) + } + + if (plot_sims == TRUE) { + DAISIE_plot_sims( + island_replicates = island_replicates, + sample_freq = sample_freq, + trait_pars = trait_pars + ) + } + return(island_replicates) +} diff --git a/R/DAISIE_sim_update_state_trait_dep.R b/R/DAISIE_sim_update_state_trait_dep.R index 6b6cf889..9ae62d23 100644 --- a/R/DAISIE_sim_update_state_trait_dep.R +++ b/R/DAISIE_sim_update_state_trait_dep.R @@ -416,93 +416,4 @@ DAISIE_sim_update_state_trait_dep <- function(timeval, return(updated_state) } -# !POTENTIALLY WRONG DUPLICATE FUNCTION! --------------------------------------- -# DAISIE_ONEcolonist <- function(time,island_spec,stt_table) -# { -# total_time <- time -# ### number of independent colonisations -# uniquecolonisation <- as.numeric(unique(island_spec[,"Colonisation time (BP)"])) -# number_colonisations <- length(uniquecolonisation) -# -# ### if there is only one independent colonisation - anagenetic and cladogenetic -# #species are classed as stac=2; immigrant classed as stac=4: -# if(number_colonisations == 1) -# { -# if(island_spec[1,"Species type"] == "I") -# { -# descendants <- list(stt_table = stt_table, -# branching_times = c(total_time,as.numeric(island_spec[1,"Colonisation time (BP)"])), -# stac = 4, -# missing_species = 0) -# } -# if(island_spec[1,"Species type"] == "A") -# { -# descendants <- list(stt_table = stt_table, -# branching_times = c(total_time,as.numeric(island_spec[1,"Colonisation time (BP)"])), -# stac = 2, -# missing_species = 0) -# } -# if(island_spec[1,"Species type"] == "C") -# { -# descendants <- list(stt_table = stt_table, -# branching_times = c(total_time,rev(sort(as.numeric(island_spec[,"branching time (BP)"])))), -# stac = 2, -# missing_species = 0) -# } -# } -# -# ### if there are two or more independent colonisations, all species are classed as stac=3 and put within same list item: -# else if(number_colonisations > 1) -# { -# descendants <- list(stt_table = stt_table, -# branching_times = NA,stac = 2,missing_species = 0, -# other_clades_same_ancestor = list()) -# ### create table with information on other clades with same ancestor, but this information is not used in DAISIE_ML -# oldest <- which(as.numeric(island_spec[,"Colonisation time (BP)"]) == max(as.numeric(island_spec[,"Colonisation time (BP)"]))) -# -# oldest_table <- island_spec[oldest,] -# if(class(oldest_table) == 'character') -# { -# oldest_table <- t(as.matrix(oldest_table)) -# } -# if(oldest_table[1,'Species type'] == 'A') -# { -# descendants$branching_times <- c(total_time, as.numeric(oldest_table[1,"Colonisation time (BP)"])) -# } else if(oldest_table[1,'Species type'] == 'C') -# { -# descendants$branching_times <- c(total_time, rev(sort(as.numeric(oldest_table[,'branching time (BP)'])))) -# } -# -# youngest_table = island_spec[-oldest,] -# if(class(youngest_table) == 'character') -# { -# youngest_table <- t(as.matrix(youngest_table)) -# } -# -# uniquecol <- as.numeric(unique(youngest_table[,"Colonisation time (BP)"])) -# -# descendants$missing_species <- length(which(youngest_table[,"Species type"]!='I')) -# for(colonisation in 1:length(uniquecol)) -# { -# descendants$other_clades_same_ancestor[[colonisation]] <- list(brts_miss = NA,species_type = NA) -# -# samecolonisation <- which(as.numeric(youngest_table[,"Colonisation time (BP)"]) == uniquecol[colonisation]) -# -# if(youngest_table[samecolonisation[1],"Species type"] == "I") -# { -# descendants$stac <- 3 -# descendants$other_clades_same_ancestor[[colonisation]]$brts_miss <- as.numeric(youngest_table[samecolonisation,"Colonisation time (BP)"]) -# descendants$other_clades_same_ancestor[[colonisation]]$species_type <- "I" -# } else if(youngest_table[samecolonisation[1],"Species type"] == "A") -# { -# descendants$other_clades_same_ancestor[[colonisation]]$brts_miss <- as.numeric(youngest_table[samecolonisation,"Colonisation time (BP)"]) -# descendants$other_clades_same_ancestor[[colonisation]]$species_type <- "A" -# } else if (youngest_table[samecolonisation[1],"Species type"] == "C") -# { -# descendants$other_clades_same_ancestor[[colonisation]]$brts_miss <- rev(sort(as.numeric(youngest_table[samecolonisation,"branching time (BP)"]))) -# descendants$other_clades_same_ancestor[[colonisation]]$species_type <- "C" -# } -# } -# } -# descendants -# } + diff --git a/R/create_pars.R b/R/create_pars.R index b2576576..19eac88e 100644 --- a/R/create_pars.R +++ b/R/create_pars.R @@ -121,20 +121,20 @@ create_trait_pars <- function(trans_rate, clado_rate2, trans_rate2, M2) { - testit::assert(is.numeric(trans_rate)) - testit::assert(is.numeric(immig_rate2)) - testit::assert(is.numeric(ext_rate2)) - testit::assert(is.numeric(ana_rate2)) - testit::assert(is.numeric(clado_rate2)) - testit::assert(is.numeric(trans_rate2)) - testit::assert(floor(M2) == M2) - testit::assert(trans_rate >= 0.0) - testit::assert(immig_rate2 >= 0.0) - testit::assert(ext_rate2 >= 0.0) - testit::assert(ana_rate2 >= 0.0) - testit::assert(clado_rate2 >= 0.0) - testit::assert(trans_rate2 >=0.0) - testit::assert(M2 >=0) + # testit::assert(is.numeric(trans_rate)) + # testit::assert(is.numeric(immig_rate2)) + # testit::assert(is.numeric(ext_rate2)) + # testit::assert(is.numeric(ana_rate2)) + # testit::assert(is.numeric(clado_rate2)) + # testit::assert(is.numeric(trans_rate2)) + # testit::assert(floor(M2) == M2) + # testit::assert(trans_rate >= 0.0) + # testit::assert(immig_rate2 >= 0.0) + # testit::assert(ext_rate2 >= 0.0) + # testit::assert(ana_rate2 >= 0.0) + # testit::assert(clado_rate2 >= 0.0) + # testit::assert(trans_rate2 >=0.0) + # testit::assert(M2 >=0) list(trans_rate = trans_rate, immig_rate2 = immig_rate2, ext_rate2 = ext_rate2, @@ -144,6 +144,41 @@ create_trait_pars <- function(trans_rate, M2 = M2) } + + +#' Create named list of trait state parameters +#' +#' @param trans_rate A numeric with the per capita transition rate with state1 +#' @param immig_rate2 A numeric with the per capita immigration rate with state2 +#' @param ext_rate2 A numeric with the per capita extinction rate with state2 +#' @param ana_rate2 A numeric with the per capita anagenesis rate with state2 +#' @param clado_rate2 A numeric with the per capita cladogenesis rate with state2 +#' @param trans_rate2 A numeric with the per capita transition rate with state2 +#' @param M2 A numeric with the number of species with trait state 2 on mainland +#' @param K2 A numeric with the carrying capacity for state 2 +#' +#' @return list of numerical values containing trait state parameters +#' @export +#' +create_trait_pars_2K <- function(trans_rate, + immig_rate2, + ext_rate2, + ana_rate2, + clado_rate2, + trans_rate2, + M2, + K2) { + list(trans_rate = trans_rate, + immig_rate2 = immig_rate2, + ext_rate2 = ext_rate2, + ana_rate2 = ana_rate2, + clado_rate2 = clado_rate2, + trans_rate2 = trans_rate2, + M2 = M2, + K2 = K2) +} + + #' Creates the list object for CS_version argument in DAISIE_ML_CS #' #' @param model the CS model to run, options are \code{1} for single rate @@ -152,6 +187,8 @@ create_trait_pars <- function(trans_rate, #' @param relaxed_par the parameter to relax (integrate over). Options are #' \code{"cladogenesis"}, \code{"extinction"}, \code{"carrying_capacity"}, #' \code{"immigration"}, or \code{"anagenesis"} +#' @param par_sd standard deviation of the parameter to relax +#' @param par_upper_bound upper bound of the parameter to relax #' @return A list of two elements #' \itemize{ #' \item{model: the CS model to run, options are \code{1} for single rate @@ -163,7 +200,9 @@ create_trait_pars <- function(trans_rate, #' } #' @export create_CS_version <- function(model = 1, - relaxed_par = NULL) { + relaxed_par = NULL, + par_sd = 0, + par_upper_bound = Inf) { if (model != 1 && model != 2 && model != 3) { stop("model must be either 1, 2 or 3") @@ -172,7 +211,9 @@ create_CS_version <- function(model = 1, stop("relaxed_par required for multi-rate model") } CS_version <- list(model = model, - relaxed_par = relaxed_par) + relaxed_par = relaxed_par, + par_sd = par_sd, + par_upper_bound = par_upper_bound) return(CS_version) } diff --git a/R/default_params_doc.R b/R/default_params_doc.R index f8381be2..e205477a 100644 --- a/R/default_params_doc.R +++ b/R/default_params_doc.R @@ -469,6 +469,12 @@ #' @param carr_cap Numeric carrying capacity #' @param immig_rate Numeric rate of immigration #' @param ana_rate Numeric rate of anagenesis +#' @param islands Island datalist or simulated data in DAISIE datalist format. +#' Can be a single island (empirical data) generated with DAISIE_dataprep or +#' DAISIEprep. Can also be simulated data generated with DAISIE_sim function. +#' @param sort_clade_sizes Default sort_clade_sizes=T outputs clade sizes +#' sorted in ascending order of number of species. sort_clade_sizes=F outputs +#' clade sizes in the same order as they appear in the input datalist. #' #' #' @return Nothing diff --git a/README.md b/README.md index e01d37ef..e4347b4c 100644 --- a/README.md +++ b/README.md @@ -86,29 +86,28 @@ For feature requests or bug-reports or other matters, please submit an [issue](h ## References -Etienne R. S., Valente, L., Phillimore, A. B., Haegeman, B., Lambert, J. W., Neves, P., Xie, S., Bilderbeek, R. J. C., & Hildenbrandt, H. (2021). DAISIE: Dynamical Assembly of Islands by Speciation, Immigration and Extinction. R package version 4.0.2. https://cran.r-project.org/package=DAISIE. https://doi.org/10.5281/zenodo.4054058 +* Etienne R. S., Valente, L., Phillimore, A. B., Haegeman, B., Lambert, J. W., Neves, P., Xie, S., Bilderbeek, R. J. C., & Hildenbrandt, H. (2021). DAISIE: Dynamical Assembly of Islands by Speciation, Immigration and Extinction. R package version 4.0.2. https://cran.r-project.org/package=DAISIE. https://doi.org/10.5281/zenodo.4054058 -Valente, L., Etienne, R.S., & Phillimore, A.B. (2014). The effects of island ontogeny on species diversity and phylogeny. Proceedings of the Royal Society B: Biological Sciences, 281(1784), 20133227–20133227. http://doi.org/10.1098/rspb.2013.3227 +* Valente, L., Etienne, R.S., & Phillimore, A.B. (2014). The effects of island ontogeny on species diversity and phylogeny. Proceedings of the Royal Society B: Biological Sciences, 281(1784), 20133227–20133227. http://doi.org/10.1098/rspb.2013.3227 -Valente, L., Phillimore, A.B., & Etienne, R.S. (2015). Equilibrium and non-equilibrium dynamics simultaneously operate in the Galápagos islands. Ecology Letters, 18(8), 844–852. http://doi.org/10.1111/ele.12461 +* Valente, L., Phillimore, A.B., & Etienne, R.S. (2015). Equilibrium and non-equilibrium dynamics simultaneously operate in the Galápagos islands. Ecology Letters, 18(8), 844–852. http://doi.org/10.1111/ele.12461 -Valente, L., Etienne, R.S., & Dávalos, L.M. (2017). Recent extinctions disturb path to equilibrium diversity in Caribbean bats. Nature Ecology & Evolution, 1(2), 0026. http://doi.org/10.1038/s41559-016-0026 +* Valente, L., Etienne, R.S., & Dávalos, L.M. (2017). Recent extinctions disturb path to equilibrium diversity in Caribbean bats. Nature Ecology & Evolution, 1(2), 0026. http://doi.org/10.1038/s41559-016-0026 -Valente, L., Illera, J.C., Havenstein, K., Pallien, T., Etienne, R.S., & Tiedemann, R. (2017). Equilibrium Bird Species Diversity in Atlantic Islands. Current Biology, 27(11), 1660-1666. https://doi.org/10.1016/j.cub.2017.04.053 +* Valente, L., Illera, J.C., Havenstein, K., Pallien, T., Etienne, R.S., & Tiedemann, R. (2017). Equilibrium Bird Species Diversity in Atlantic Islands. Current Biology, 27(11), 1660-1666. https://doi.org/10.1016/j.cub.2017.04.053 -Valente, L., Phillimore, A.B., & Etienne, R.S. (2018). Using molecular phylogenies in island biogeography: It’s about time. Ecography, 1–3. http://doi.org/10.1111/ecog.03503 +* Valente, L., Phillimore, A.B., & Etienne, R.S. (2018). Using molecular phylogenies in island biogeography: It’s about time. Ecography, 1–3. http://doi.org/10.1111/ecog.03503 -Valente, L., Etienne, R.S., & Garcia-R., J.C. (2019). Deep Macroevolutionary Impact of Humans on New Zealand’s Unique Avifauna. Current Biology 29 (15): 2563-2569.e4. https://doi.org/10.1016/j.cub.2019.06.058 +* Valente, L., Etienne, R.S., & Garcia-R., J.C. (2019). Deep Macroevolutionary Impact of Humans on New Zealand’s Unique Avifauna. Current Biology 29 (15): 2563-2569.e4. https://doi.org/10.1016/j.cub.2019.06.058 -Valente, L., Phillimore, A.B., Melo, M., Warren, B.H., Clegg, S.M., Havenstein, K., Tiedemann, R., Illera, J.C., Thebaud, C., Aschenbach, T. & Etienne, R.S. (2020). A Simple Dynamic Model Explains the Diversity of Island Birds Worldwide. Nature 579 (7797): 92–96. https://doi.org/10.1038/s41586-020-2022-5 +* Valente, L., Phillimore, A.B., Melo, M., Warren, B.H., Clegg, S.M., Havenstein, K., Tiedemann, R., Illera, J.C., Thebaud, C., Aschenbach, T. & Etienne, R.S. (2020). A Simple Dynamic Model Explains the Diversity of Island Birds Worldwide. Nature 579 (7797): 92–96. https://doi.org/10.1038/s41586-020-2022-5 -Hauffe, T., Delicado, D., Etienne, R.S., & Valente, L. (2020). Lake expansion elevates equilibrium diversity via increasing colonization. Journal of Biogeography 47: 1849–1860. https://doi.org/10.1111/jbi.13914 +* Hauffe, T., Delicado, D., Etienne, R.S., & Valente, L. (2020). Lake expansion elevates equilibrium diversity via increasing colonization. Journal of Biogeography 47: 1849–1860. https://doi.org/10.1111/jbi.13914 -Valente, L., Kristensen, N., Phillimore, A. B., & Etienne, R. S. (2021). Report of programming bugs in the DAISIE R package: consequences and correction. https://doi.org/10.32942/osf.io/w5ntf +* Valente, L., Kristensen, N., Phillimore, A. B., & Etienne, R. S. (2021). Report of programming bugs in the DAISIE R package: consequences and correction. https://doi.org/10.32942/osf.io/w5ntf +* Santos Neves, P.\*, Lambert, J. W.\*, Valente, L., & Etienne, R. S. (2022). The robustness of a simple dynamic model of island biodiversity to geological and sea-level change. Journal of Biogeography. https://doi.org/10.1111/jbi.14519 -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. https://doi.org/10.1101/2021.07.26.453064 +* Lambert, J. W., Santos Neves, P., Bilderbeek, R. L. C., Valente, L., Etienne, R. S. (2022). The effect of mainland dynamics on data and parameter estimates in island biogeography. bioRxiv. https://doi.org/10.1101/2022.01.13.476210 -Lambert, J. W., Santos Neves, P., Bilderbeek, R. L. C., Valente, L., Etienne, R. S. (2022). The effect of mainland dynamics on data and parameter estimates in island biogeography. bioRxiv. https://doi.org/10.1101/2022.01.13.476210 - -Xie, S., Valente, L., Etienne, R. S. (2022). A simple island biodiversity model is robust to trait dependence in diversification and colonization rates. biRrxiv. https://doi.org/10.1101/2022.01.01.474685 +* Xie, S., Valente, L., Etienne, R. S. (2023). Can we ignore trait-dependent colonization and diversification in island biogeography? Evolution. https://doi.org/10.1093/evolut/qpad006. diff --git a/man/DAISIE-package.Rd b/man/DAISIE-package.Rd index 9260d5d4..464c1708 100644 --- a/man/DAISIE-package.Rd +++ b/man/DAISIE-package.Rd @@ -23,9 +23,9 @@ Cladogenesis and immigration rates can be dependent on diversity. \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 \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. \doi{10.1101/2021.07.26.453064}.\cr +\item Santos Neves, P., Lambert, J. W., Valente, L., & Etienne, R. S. (2022). The robustness of a simple dynamic model of island biodiversity to geological and sea-level change. Journal of Biogeography. \doi{10.1111/jbi.14519}.\cr \item Lambert, J. W., Santos Neves, P., Bilderbeek, R. L. C., Valente, L., Etienne, R. S. (2022). The effect of mainland dynamics on data and parameter estimates in island biogeography. bioRxiv. \doi{10.1101/2022.01.13.476210}.\cr -\item Xie, S., Valente, L., Etienne, R. S. (2022). A simple island biodiversity model is robust to trait dependence in diversification and colonization rates. biRrxiv. \doi{10.1101/2022.01.01.474685}.\cr +\item Xie, S., Valente, L., Etienne, R. S. (2023). Can we ignore trait-dependent colonization and diversification in island biogeography? Evolution. \doi{10.1093/evolut/qpad006}.\cr } } \seealso{ @@ -45,10 +45,10 @@ Authors: \item Albert B. Phillimore (\href{https://orcid.org/0000-0002-6553-1553}{ORCID}) \item Bart Haegeman (\href{https://orcid.org/0000-0003-2325-4727}{ORCID}) \item Joshua W. Lambert \email{j.w.l.lambert@rug.nl} (\href{https://orcid.org/0000-0001-5218-3046}{ORCID}) - \item Pedro Neves \email{p.m.santos.neves@rug.nl} (\href{https://orcid.org/0000-0003-2561-4677}{ORCID}) + \item Pedro Santos Neves \email{p.m.santos.neves@rug.nl} (\href{https://orcid.org/0000-0003-2561-4677}{ORCID}) \item Shu Xie \email{s.xie@rug.nl} (\href{https://orcid.org/0000-0001-9594-946X}{ORCID}) \item Richèl J.C. Bilderbeek \email{richel@richelbilderbeek.nl} (\href{https://orcid.org/0000-0003-1107-7049}{ORCID}) - \item Hanno Hildenbrandt \email{h.hildenbrandt@rug.nl} + \item Hanno Hildenbrandt \email{h.hildenbrandt@rug.nl} (\href{https://orcid.org/0000-0002-6784-1037}{ORCID}) } Other contributors: diff --git a/man/DAISIE_ML4.Rd b/man/DAISIE_ML4.Rd index 78c068d8..b30ae271 100644 --- a/man/DAISIE_ML4.Rd +++ b/man/DAISIE_ML4.Rd @@ -18,7 +18,8 @@ DAISIE_ML4( maxiter = 1000 * round((1.25)^length(idparsopt)), methode = "lsodes", optimmethod = "subplex", - CS_version = create_CS_version(model = 2, relaxed_par = "cladogenesis"), + CS_version = create_CS_version(model = 2, relaxed_par = "cladogenesis", par_sd = 0, + par_upper_bound = Inf), verbose = 0, tolint = c(1e-16, 1e-10), island_ontogeny = NA, diff --git a/man/DAISIE_count_species.Rd b/man/DAISIE_count_species.Rd new file mode 100644 index 00000000..c522d1eb --- /dev/null +++ b/man/DAISIE_count_species.Rd @@ -0,0 +1,59 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/DAISIE_count_species.R +\name{DAISIE_count_species} +\alias{DAISIE_count_species} +\title{Count number of species in DAISIE datalist or simulated data.} +\usage{ +DAISIE_count_species(islands, sort_clade_sizes = TRUE) +} +\arguments{ +\item{islands}{Island datalist or simulated data in DAISIE datalist format. +Can be a single island (empirical data) generated with DAISIE_dataprep or +DAISIEprep. Can also be simulated data generated with DAISIE_sim function.} + +\item{sort_clade_sizes}{Default sort_clade_sizes=T outputs clade sizes +sorted in ascending order of number of species. sort_clade_sizes=F outputs +clade sizes in the same order as they appear in the input datalist.} +} +\value{ +The output is a list containing the following items: +\item{clade_sizes_sorted}{ List showing the total number of species in each +island clade (including missing species). Each item [[i]] on the list +gives the sizes of all clades for a single island. If option +sort_clade_sizes = T, +the clade sizes for are +sorted by increasing number of species. If option sort_clade_sizes = F +the clade sizes are given in the same order as in the input datalist.} +\item{size_largest_clade}{ The total number of species in the largest +island clade +for each island.} +\item{mean_clade_size}{ Mean clade size (average of all island clades)} +\item{number_colonisations}{ The total number of colonisations (clades) on +each island.} +\item{total_number_species}{ The total number of species on each island. These +are the extant species at present, including missing species; in case of +simulations, this is the number of species present on the island at the +end of the +simulation.} +} +\description{ +Calculates various island diversity metrics from island datasets. +} +\examples{ +# Run function with clade sizes in the order they appear in the input data +data("NewZealand_birds_datalist") +species_count <- DAISIE_count_species(NewZealand_birds_datalist) + +# Run function with clade sizes in ascending order +species_count_sorted <- DAISIE_count_species( + NewZealand_birds_datalist, + sort_clade_sizes = TRUE +) +} +\seealso{ +\code{\link{DAISIE_dataprep}}, +\code{\link{DAISIE_plot_island}} +} +\author{ +Luis Valente +} diff --git a/man/DAISIE_sim_relaxed_rate.Rd b/man/DAISIE_sim_relaxed_rate.Rd index 6158fd41..fc62a68b 100644 --- a/man/DAISIE_sim_relaxed_rate.Rd +++ b/man/DAISIE_sim_relaxed_rate.Rd @@ -183,8 +183,8 @@ ext_rate <- 0.2 carr_cap <- Inf immig_rate <- 0.005 ana_rate <- 1 -sd <- 1 -sim_pars <- c(clado_rate, ext_rate, carr_cap, immig_rate, ana_rate, sd) +par_sd <- 1 +sim_pars <- c(clado_rate, ext_rate, carr_cap, immig_rate, ana_rate, par_sd) set.seed(1) island_replicates <- DAISIE_sim_relaxed_rate( time = 1, diff --git a/man/DAISIE_sim_trait_dep.Rd b/man/DAISIE_sim_trait_dep.Rd index fada74dc..6e082fe5 100644 --- a/man/DAISIE_sim_trait_dep.Rd +++ b/man/DAISIE_sim_trait_dep.Rd @@ -2,7 +2,7 @@ % Please edit documentation in R/DAISIE_sim_trait_dep.R \name{DAISIE_sim_trait_dep} \alias{DAISIE_sim_trait_dep} -\title{Simulate islands with given parameters.} +\title{Simulate islands with given trait-dependent parameters.} \usage{ DAISIE_sim_trait_dep( time, @@ -10,7 +10,6 @@ DAISIE_sim_trait_dep( pars, replicates, divdepmodel = "CS", - num_guilds = NULL, sample_freq = 25, plot_sims = TRUE, island_ontogeny = "const", @@ -71,10 +70,6 @@ carrying capacity, where diversity-dependence operates within and among clades. Option divdepmodel = 'GW' runs a model with diversity-dependence operates within a guild.} -\item{num_guilds}{The number of guilds on the mainland. The number of -mainland species is divided by the number of guilds when -\code{divdepmodel = "GW"}} - \item{sample_freq}{Numeric specifing the number of units times should be divided by for plotting purposes. Larger values will lead to plots with higher resolution, but will also run slower.} @@ -156,52 +151,39 @@ means also intermediate progress during loglikelihood computation is shown.} \item{...}{Any arguments to pass on to plotting functions.} } \value{ -Each simulated dataset is an element of the list, which can be -called using [[x]]. For example if the object is called island_replicates, -the last replicates is a list in itself. The first (e.g. -\code{island_replicates[[x]][[1]]}) element of that list has the following -components: \cr \code{$island_age} - the island age \cr Then, depending on -whether a distinction between types is made, we have:\cr \code{$not_present} -- the number of mainland lineages that are not present on the island \cr -or:\cr \code{$not_present_type1} - the number of mainland lineages of type 1 -that are not present on the island \cr \code{$not_present_type2} - the -number of mainland lineages of type 2 that are not present on the island \cr -\code{$stt_all} - STT table for all species on the island (nI - number of -non-endemic species; nA - number of anagenetic species, nC - number of -cladogenetic species, present - number of independent colonisations present -)\cr \code{$stt_stt_type1} - STT table for type 1 species on the island - -only if 2 types of species were simulated (nI - number of non-endemic -species; nA - number of anagenetic species, nC - number of cladogenetic -species, present - number of independent colonisations present )\cr -\code{$stt_stt_type2} - STT table for type 2 species on the island - only if -2 types of species were simulated (nI - number of non-endemic species; nA - -number of anagenetic species, nC - number of cladogenetic species, present - -number of independent colonisations present )\cr \code{$brts_table} - Only -for simulations under 'IW'. Table containing information on order of events -in the data, for use in maximum likelihood optimization.)\cr - -The subsequent elements of the list each contain information on a single -colonist lineage on the island and has 4 components:\cr -\code{$branching_times} - island age and stem age of the population/species -in the case of Non-endemic, Non-endemic_MaxAge and Endemic anagenetic -species. For cladogenetic species these should be island age and branching -times of the radiation including the stem age of the radiation.\cr -\code{$stac} - the status of the colonist \cr * Non_endemic_MaxAge: 1 \cr * -ndemic: 2 \cr * Endemic&Non_Endemic: 3 \cr * Non_endemic: 4 \cr -\code{$missing_species} - number of island species that were not sampled for -particular clade (only applicable for endemic clades) \cr \code{$type_1or2} -- whether the colonist belongs to type 1 or type 2 \cr +A list. The highest level of the least corresponds to each individual +replciate. The first element of each replicate is composed of island +information containing: +\itemize{ + \item{\code{$island_age}: A numeric with the island age.} + \item{\code{$not_present}: A numeric with the number of mainland lineages + that are not present on the island.} + \item{\code{$stt_all}: STT table for all species on the island + (nI - number of non-endemic species; nA - number of anagenetic species, + nC - number of cladogenetic species, present - number of independent + colonisations present)} + \item{\code{$brts_table}: Only for simulations under \code{"IW"}. Table +containing information on order of events in the data, for use in maximum +likelihood optimization.).} +} +The subsequent elements of the list pertaining to each replcate contain +information on a single colonist lineage on the island and have 4 components: +\itemize{ + \item{\code{$branching_times}: island age and stem age of the + population/species in the case of Non-endemic, Non-endemic_MaxAge and + Endemic anagenetic species. For cladogenetic species these should + be island age and branching times of the radiation including the + stem age of the radiation.} + \item{\code{$stac}: An integer ranging from 1 to 4 + indicating the status of the colonist:} +\item{\code{$missing_species}: number of island species that were +not sampled for particular clade (only applicable for endemic clades)} +\item{\code{$type_1or2}: whether the colonist belongs to type 1 or type 2} +} } \description{ This function simulates islands with given cladogenesis, -extinction, Kprime, immigration and anagenesis parameters. If a single -parameter set is provided (5 parameters) it simulates islands where all -species have the same macro-evolutionary process. If two paramater sets -(10 parameters) are provided, it simulates islands where two different -macro-evolutionary processes operate, one applying to type 1 species and -other to type 2 species. If two parameter sets and a time shift (11 -parameters) are provided, it simulates islands where at the given time -a shift between the parameter sets will occur. +extinction, K, immigration and anagenesis parameters for binary states. Returns R list object that contains the simulated islands } diff --git a/man/DAISIE_sim_trait_dep_2K.Rd b/man/DAISIE_sim_trait_dep_2K.Rd new file mode 100644 index 00000000..414661a5 --- /dev/null +++ b/man/DAISIE_sim_trait_dep_2K.Rd @@ -0,0 +1,205 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/DAISIE_sim_trait_dep_2K.R +\name{DAISIE_sim_trait_dep_2K} +\alias{DAISIE_sim_trait_dep_2K} +\title{Simulate islands with given trait-dependent parameters.} +\usage{ +DAISIE_sim_trait_dep_2K( + time, + M, + pars, + replicates, + divdepmodel = "CS", + sample_freq = 25, + plot_sims = TRUE, + island_ontogeny = "const", + sea_level = "const", + hyper_pars = create_hyper_pars(d = 0, x = 0), + area_pars = DAISIE::create_area_pars(max_area = 1, current_area = 1, + proportional_peak_t = 0, total_island_age = 0, sea_level_amplitude = 0, + sea_level_frequency = 0, island_gradient_angle = 0), + extcutoff = 1000, + cond = 0, + verbose = TRUE, + trait_pars = NULL, + ... +) +} +\arguments{ +\item{time}{Numeric defining the length of the simulation in time units. +For example, if an island is known to be 4 million years old, setting +time = 4 will simulate the entire life span of the island; setting time = 2 +will stop the simulation at the mid-life of the island.} + +\item{M}{Numeric defining the size of mainland pool, i.e. the number of +species that can potentially colonize the island.} + +\item{pars}{A numeric vector containing the model parameters: +\itemize{ + \item{\code{pars[1]}: lambda^c (cladogenesis rate)} + \item{\code{pars[2]}: mu (extinction rate)} + \item{\code{pars[3]}: K (carrying capacity), set K=Inf for diversity + independence.} + \item{\code{pars[4]}: gamma (immigration rate)} + \item{\code{pars[5]}: lambda^a (anagenesis rate)} + \item{\code{pars[6]}: lambda^c (cladogenesis rate) for either type 2 species + or rate set 2 in rate shift model} + \item{\code{pars[7]}: mu (extinction rate) for either type 2 species or rate + set 2 in rate shift model} + \item{\code{pars[8]}: K (carrying capacity) for either type 2 species or rate + set 2 in rate shift model, set K=Inf for diversity independence.} + \item{\code{pars[9]}: gamma (immigration rate) for either type 2 species + or rate set 2 in rate shift model} + \item{\code{pars[10]}: lambda^a (anagenesis rate) for either type 2 + species or rate set 2 in rate shift model} +} +Elements 6:10 are required only when type 2 species are included +or in the rate shift model. For \code{\link{DAISIE_sim_relaxed_rate}()} +\code{pars[6]} is the standard deviation of the gamma distribution for the +relaxed parameter and the parameter chosen by the \code{relaxed_par} +argument is the mean of the gamma distribution for the relaxed parameter.} + +\item{replicates}{Integer specifying number of island replicates to be +simulated.} + +\item{divdepmodel}{Option divdepmodel = 'CS' runs a model with clade-specific +carrying capacity, where diversity-dependence operates only within single +clades, i.e. only among species originating from the same mainland +colonist. Option divdepmodel = 'IW' runs a model with island-wide +carrying capacity, where diversity-dependence operates within and among +clades. Option divdepmodel = 'GW' runs a model with diversity-dependence +operates within a guild.} + +\item{sample_freq}{Numeric specifing the number of units times should be +divided by for plotting purposes. Larger values will lead to plots with +higher resolution, but will also run slower.} + +\item{plot_sims}{\code{Default = TRUE} plots species-through-time (STT) +plots. It detects how many types of species are present. If only one type +of species is present, STT is plotted for all species. If two types are +present, three plots are produced: STT for all, STT for type 1 and STT for +type 2.} + +\item{island_ontogeny}{In \code{\link{DAISIE_sim_time_dep}()}, +\code{\link{DAISIE_ML_CS}} and plotting a string describing the type of +island ontogeny. Can be \code{"const"}, \code{"beta"} for a beta function +describing area through time. String checked by +\code{\link{is_island_ontogeny_input}()}. \cr In all other functions a +numeric describing the type of island ontogeny. Can be \code{0} for +constant, \code{1} for a beta function describing area through time. In ML +functions \code{island_ontogeny = NA} assumes constant ontogeny. Time +dependent estimation is not yet available as development is still ongoing. +Will return an error if called in that case.} + +\item{sea_level}{In \code{\link{DAISIE_sim_time_dep}()} and plotting a +string describing the type of sea level. Can be \code{"const"} or +\code{"sine"} for a sine function describing area through time. String +checked by \code{\link{is_sea_level_input}()}. +\cr In all other functions a numeric describing the type of sea level. Can +be \code{0} for constant, \code{1} for a sine function describing area +through time.} + +\item{hyper_pars}{A named list of numeric hyperparameters for the rate +calculations as returned by \code{\link{create_hyper_pars}()}: +\itemize{ + \item{[1]: is d the scaling parameter for exponent for calculating + cladogenesis rate} + \item{[2]: is x the exponent for calculating extinction rate} +}} + +\item{area_pars}{A named list containing area and sea level parameters as +created by \code{\link{create_area_pars}()}: +\itemize{ + \item{[1]: maximum area} + \item{[2]: current area} + \item{[3]: value from 0 to 1 indicating where in the island's history the + peak area is achieved} + \item{[4]: total island age} + \item{[5]: amplitude of area fluctuation from sea level} + \item{[6]: frequency of sine wave of area change from sea level} + \item{[7]: angle of the slope of the island} +}} + +\item{extcutoff}{A numeric with the cutoff for the the maximum extinction +rate preventing it from being too large and slowing down simulation.} + +\item{cond}{cond = 0 : conditioning on island age \cr cond = 1 : +conditioning on island age and non-extinction of the island biota \cr. +cond > 1 : conditioning on island age and having at least cond colonizations +on the island. This last option is not yet available for the IW model \cr} + +\item{verbose}{In simulation and dataprep functions a logical, +\code{Default = TRUE} gives intermediate output should be printed. +For ML functions a numeric determining if intermediate output should be +printed, \code{Default = 0} does not print, \code{verbose = 1} prints +intermediate output of the parameters and loglikelihood, \code{verbose = 2} +means also intermediate progress during loglikelihood computation is shown.} + +\item{trait_pars}{A named list containing diversification rates considering + two trait states created by \code{\link{create_trait_pars}}: +\itemize{ + \item{[1]:A numeric with the per capita transition rate with state1} + \item{[2]:A numeric with the per capita immigration rate with state2} + \item{[3]:A numeric with the per capita extinction rate with state2} + \item{[4]:A numeric with the per capita anagenesis rate with state2} + \item{[5]:A numeric with the per capita cladogenesis rate with state2} + \item{[6]:A numeric with the per capita transition rate with state2} + \item{[7]:A numeric with the number of species with trait state 2 on + mainland} +}} + +\item{...}{Any arguments to pass on to plotting functions.} +} +\value{ +A list. The highest level of the least corresponds to each individual +replciate. The first element of each replicate is composed of island +information containing: +\itemize{ + \item{\code{$island_age}: A numeric with the island age.} + \item{\code{$not_present}: A numeric with the number of mainland lineages + that are not present on the island.} + \item{\code{$stt_all}: STT table for all species on the island + (nI - number of non-endemic species; nA - number of anagenetic species, + nC - number of cladogenetic species, present - number of independent + colonisations present)} + \item{\code{$brts_table}: Only for simulations under \code{"IW"}. Table +containing information on order of events in the data, for use in maximum +likelihood optimization.).} +} +The subsequent elements of the list pertaining to each replcate contain +information on a single colonist lineage on the island and have 4 components: +\itemize{ + \item{\code{$branching_times}: island age and stem age of the + population/species in the case of Non-endemic, Non-endemic_MaxAge and + Endemic anagenetic species. For cladogenetic species these should + be island age and branching times of the radiation including the + stem age of the radiation.} + \item{\code{$stac}: An integer ranging from 1 to 4 + indicating the status of the colonist:} +\item{\code{$missing_species}: number of island species that were +not sampled for particular clade (only applicable for endemic clades)} +\item{\code{$type_1or2}: whether the colonist belongs to type 1 or type 2} +} +} +\description{ +This function simulates islands with given cladogenesis, +extinction, K, immigration and anagenesis parameters. In this version, +rates and K are both trait-dependent. + +Returns R list object that contains the simulated islands +} +\references{ +Valente, L.M., A.B. Phillimore and R.S. Etienne (2015). +Equilibrium and non-equilibrium dynamics simultaneously operate in the +Galapagos islands. Ecology Letters 18: 844-852. +Hauffe, T., D. Delicado, R.S. Etienne and L. Valente (submitted). +Lake expansion increases equilibrium diversity via the target effect of +island biogeography. +} +\seealso{ +\code{\link{DAISIE_format_CS}} \code{\link{DAISIE_plot_sims}} +} +\author{ +Luis Valente and Albert Phillimore +} +\keyword{models} diff --git a/man/create_CS_version.Rd b/man/create_CS_version.Rd index df1864ff..0409316d 100644 --- a/man/create_CS_version.Rd +++ b/man/create_CS_version.Rd @@ -4,7 +4,12 @@ \alias{create_CS_version} \title{Creates the list object for CS_version argument in DAISIE_ML_CS} \usage{ -create_CS_version(model = 1, relaxed_par = NULL) +create_CS_version( + model = 1, + relaxed_par = NULL, + par_sd = 0, + par_upper_bound = Inf +) } \arguments{ \item{model}{the CS model to run, options are \code{1} for single rate @@ -14,6 +19,10 @@ model} \item{relaxed_par}{the parameter to relax (integrate over). Options are \code{"cladogenesis"}, \code{"extinction"}, \code{"carrying_capacity"}, \code{"immigration"}, or \code{"anagenesis"}} + +\item{par_sd}{standard deviation of the parameter to relax} + +\item{par_upper_bound}{upper bound of the parameter to relax} } \value{ A list of two elements diff --git a/man/create_trait_pars_2K.Rd b/man/create_trait_pars_2K.Rd new file mode 100644 index 00000000..da15473f --- /dev/null +++ b/man/create_trait_pars_2K.Rd @@ -0,0 +1,40 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/create_pars.R +\name{create_trait_pars_2K} +\alias{create_trait_pars_2K} +\title{Create named list of trait state parameters} +\usage{ +create_trait_pars_2K( + trans_rate, + immig_rate2, + ext_rate2, + ana_rate2, + clado_rate2, + trans_rate2, + M2, + K2 +) +} +\arguments{ +\item{trans_rate}{A numeric with the per capita transition rate with state1} + +\item{immig_rate2}{A numeric with the per capita immigration rate with state2} + +\item{ext_rate2}{A numeric with the per capita extinction rate with state2} + +\item{ana_rate2}{A numeric with the per capita anagenesis rate with state2} + +\item{clado_rate2}{A numeric with the per capita cladogenesis rate with state2} + +\item{trans_rate2}{A numeric with the per capita transition rate with state2} + +\item{M2}{A numeric with the number of species with trait state 2 on mainland} + +\item{K2}{A numeric with the carrying capacity for state 2} +} +\value{ +list of numerical values containing trait state parameters +} +\description{ +Create named list of trait state parameters +} diff --git a/man/cum_rho.Rd b/man/cum_rho.Rd new file mode 100644 index 00000000..10257ffc --- /dev/null +++ b/man/cum_rho.Rd @@ -0,0 +1,19 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/DAISIE_loglik_integrate.R +\name{cum_rho} +\alias{cum_rho} +\title{Cumulative Gamma distribution parameterised with mean and standard deviation} +\usage{ +cum_rho(par_upper_bound, DAISIE_dist_pars) +} +\arguments{ +\item{DAISIE_dist_pars}{A numeric vector of two elements, first is the mean +and second the standard deviation of the distribution.} +} +\value{ +Numeric +} +\description{ +Cumulative Gamma distribution parameterised with mean and standard deviation +} +\keyword{internal} diff --git a/man/default_params_doc.Rd b/man/default_params_doc.Rd index 1cd4d25c..3c74d4a7 100644 --- a/man/default_params_doc.Rd +++ b/man/default_params_doc.Rd @@ -726,6 +726,14 @@ relaxed-rate model} \item{immig_rate}{Numeric rate of immigration} \item{ana_rate}{Numeric rate of anagenesis} + +\item{islands}{Island datalist or simulated data in DAISIE datalist format. +Can be a single island (empirical data) generated with DAISIE_dataprep or +DAISIEprep. Can also be simulated data generated with DAISIE_sim function.} + +\item{sort_clade_sizes}{Default sort_clade_sizes=T outputs clade sizes +sorted in ascending order of number of species. sort_clade_sizes=F outputs +clade sizes in the same order as they appear in the input datalist.} } \value{ Nothing diff --git a/man/integral_peak.Rd b/man/integral_peak.Rd index 86f33bcb..0a837384 100644 --- a/man/integral_peak.Rd +++ b/man/integral_peak.Rd @@ -21,7 +21,8 @@ integral_peak( verbose, pick, par_mean, - par_sd + par_sd, + par_upper_bound ) } \arguments{ diff --git a/src/DAISIE_CS.cpp b/src/DAISIE_CS.cpp index 902ffdeb..455d2574 100644 --- a/src/DAISIE_CS.cpp +++ b/src/DAISIE_CS.cpp @@ -7,6 +7,8 @@ #include "DAISIE_odeint.h" +using namespace daisie_odeint::jacobian_policy; + namespace { @@ -21,29 +23,6 @@ namespace { static constexpr double default_abm_factor = 0.0001; static double abm_factor = default_abm_factor; - // - class padded_vector_view - { - public: - padded_vector_view(int pad, const double* data, int n) : - data_(data), pad_(pad), n_(n) - { - } - - // return 0.0 for indices 'i' outside [pad, pad + n) - double operator[](int i) const - { - const auto ii = i - pad_; - return (ii >= 0 && i < n_) ? *(data_ + ii) : 0.0; - } - - private: - const double* data_ = nullptr; // this->operator[pad_] == *data_ - int pad_ = 0; - int n_ = 0; - }; - - // common parameter struct param_t { @@ -62,32 +41,20 @@ namespace { class cpp_daisie_cs_runmod { public: + using jacobian = const_from_linear_rhs; + cpp_daisie_cs_runmod(param_t&& par) : p_(std::move(par)) { } // odeint interface - void operator()(const state_type& x, state_type& dx, double) const + void operator()(const state_type& x, state_type& dx, double t) const { if (++p_.steps > max_cs_steps) throw std::runtime_error("cpp_daisie_cs_runmod: too many steps"); - // xx1 = c(0,0,x[1:lx],0) - // xx2 = c(0,0,x[(lx + 1):(2 * lx)],0) - // xx3 = x[2 * lx + 1] - // using padded views instead of vectors: - const auto xx1 = padded_vector_view(2, x.data().begin(), p_.lx); - const auto xx2 = padded_vector_view(2, x.data().begin() + p_.lx, p_.lx); - const auto xx3 = x[2 * p_.lx]; - - // DO I = 1, N + 4 + 2 * kk - // laavec(I) = P(I) - // lacvec(I) = P(I + N + 4 + 2 * kk) - // muvec(I) = P(I + 2 * (N + 4 + 2 * kk)) - // gamvec(I) = P(I + 3 * (N + 4 + 2 * kk)) - // nn(I) = P(I + 4 * (N + 4 + 2 * kk)) - // ENDDO - // using views instead of vectors: + const auto xx1 = padded_vector_view<2>(x.data().begin(), p_.lx); + const auto xx2 = padded_vector_view<2>(x.data().begin() + p_.lx, p_.lx); const auto chunk = p_.lx + 4 + 2 * p_.kk; const auto laavec = p_.P.data().begin(); const auto lacvec = p_.P.data().begin() + chunk; @@ -95,78 +62,41 @@ namespace { const auto gamvec = p_.P.data().begin() + 3 * chunk; const auto nn = p_.P.data().begin() + 4 * chunk; - // DO I = 3, N + 2 - // il1(I - 2) = I + kk - 1 - // il2(I - 2) = I + kk + 1 - // il3in3(I - 2) = I + kk - // il4(I - 2) = I + kk - 2 - // in1(I - 2) = I + 2 * kk - 1 - // in2ix2(I - 2) = I + 1 - // ix1(I - 2) = I - 1 - // ix3(I - 2) = I - // ix4(I - 2) = I - 2 - // ENDDO // using offsets into our views instead of vectors: - const int il1 = 2 + p_.kk - 1; - const int il2 = 2 + p_.kk + 1; - const int il3in3 = 2 + p_.kk; - const int il4 = 2 + p_.kk - 2; - const int in1 = 2 + 2 * p_.kk - 1; - const int in2ix2 = 2 + 1; - const int ix1 = 2 - 1; - const int ix3 = 2; - const int ix4 = 2 - 2; - - // DO I = 1, N - // FF1 = laavec(il1(I) + 1) * xx2(ix1(I)) - // FF1 = FF1 + lacvec(il4(I) + 1) * xx2(ix4(I)) - // FF1 = FF1 + muvec(il2(I) + 1) * xx2(ix3(I)) - // FF1 = FF1 + lacvec(il1(I)) * nn(in1(I)) * xx1(ix1(I)) - // FF1 = FF1 + muvec(il2(I)) * nn(in2ix2(I)) * xx1(in2ix2(I)) - // FFF = (muvec(il3in3(I)) + lacvec(il3in3(I))) - // FF1 = FF1 - FFF * nn(il3in3(I)) * xx1(ix3(I)) - // FF1 = FF1 - gamvec(il3in3(I)) * xx1(ix3(I)) - // dConc(I) = FF1 - // FF1 = gamvec(il3in3(I)) * xx1(ix3(I)) - // FF1 = FF1 + lacvec(il1(I) + 1) * nn(in1(I)) * xx2(ix1(I)) - // FF1 = FF1 + muvec(il2(I)+1) * nn(in2ix2(I)) * xx2(in2ix2(I)) - // FFF = muvec(il3in3(I) + 1) + lacvec(il3in3(I) + 1) - // FF1 = FF1 - FFF * nn(il3in3(I) + 1) * xx2(ix3(I)) - // FF1 = FF1 - laavec(il3in3(I) + 1) * xx2(ix3(I)) - // dConc(N + I) = FF1 - // ENDDO - // using views into output vector: + const int nil2lx = 2; + const int il1 = nil2lx + p_.kk - 1; + const int il2 = nil2lx + p_.kk + 1; + const int il3 = nil2lx + p_.kk; + const int il4 = nil2lx + p_.kk - 2; + + const int in1 = nil2lx + 2 * p_.kk - 1; + const int in2 = nil2lx + 1; + const int in3 = nil2lx + p_.kk; + + const int ix1 = nil2lx - 1; + const int ix2 = nil2lx + 1; + const int ix3 = nil2lx; + const int ix4 = nil2lx - 2; + auto dx1 = dx.data().begin(); auto dx2 = dx1 + p_.lx; + auto dx3 = dx2 + p_.lx; + for (int i = 0; i < p_.lx; ++i) { dx1[i] = laavec[il1 + i + 1] * xx2[ix1 + i] + lacvec[il4 + i + 1] * xx2[ix4 + i] + muvec[il2 + i + 1] * xx2[ix3 + i] + lacvec[il1 + i] * nn[in1 + i] * xx1[ix1 + i] - + muvec[il2 + i] * nn[in2ix2 + i] * xx1[in2ix2 + i] - - (muvec[il3in3 + i] + lacvec[il3in3 + i]) * nn[il3in3 + i] * xx1[ix3 + i] - - gamvec[il3in3 + i] * xx1[ix3 + i]; - dx2[i] = gamvec[il3in3 + i] * xx1[ix3 + i] + + muvec[il2 + i] * nn[in2 + i] * xx1[ix2 + i] + - (muvec[il3 + i] + lacvec[il3 + i]) * nn[in3 + i] * xx1[ix3 + i] + - gamvec[il3 + i] * xx1[ix3 + i]; + dx2[i] = gamvec[il3 + i] * xx1[ix3 + i] + lacvec[il1 + i + 1] * nn[in1 + i] * xx2[ix1 + i] - + muvec[il2 + i + 1] * nn[in2ix2 + i] * xx2[in2ix2 + i] - - (muvec[il3in3 + 1 + i] + lacvec[il3in3 + 1 + i]) * nn[il3in3 + i + 1] * xx2[ix3 + i] - - laavec[il3in3 + i] * xx2[ix3 + i]; + + muvec[il2 + i + 1] * nn[in2 + i] * xx2[ix2 + i] + - (muvec[il3 + 1 + i] + lacvec[il3 + 1 + i]) * nn[in3 + i + 1] * xx2[ix3 + i] + - laavec[il3 + i] * xx2[ix3 + i]; } - - // IF(kk .EQ. 1) THEN - // dConc(1) = dConc(1) + laavec(il3in3(1)) * xx3 - // dConc(2) = dConc(2) + 2 * lacvec(il3in3(1)) * xx3 - // ENDIF - // if (1 == p_.kk) { - // dx1[0] += laavec[il3in3] * xx3; - // dx2[1] += 2.0 * lacvec[il3in3] * xx3; - //} - - // FFF = laavec(il3in3(1)) + lacvec(il3in3(1)) - // FFF = FFF + gamvec(il3in3(1)) + muvec(il3in3(1)) - // dConc(2 * N + 1) = -1 * FFF * xx3 - auto dx3 = dx2 + p_.lx; - dx3[0] = -(laavec[il3in3] + lacvec[il3in3] + gamvec[il3in3] + muvec[il3in3]) * xx3; + dx3[0] = 0.0; } private: @@ -176,8 +106,10 @@ namespace { class cpp_daisie_cs_runmod_1 { public: - cpp_daisie_cs_runmod_1(param_t&& p) : - p_(p) + using jacobian = const_from_linear_rhs; + + cpp_daisie_cs_runmod_1(param_t&& par) : + p_(std::move(par)) { } @@ -186,24 +118,11 @@ namespace { { if (++p_.steps > max_cs_steps) throw std::runtime_error("cpp_daisie_cs_runmod_1: too many steps"); - // xx1 = c(0,0,x[1:lx],0) - // xx2 = c(0,0,x[(lx + 1):(2 * lx)],0) - // xx3 = c(0,0,x[(2 * lx + 1):(3 * lx)],0) - // xx4 <- c(0,0,x[(3 * lx + 1):(4 * lx)],0) - // using padded views instead of vectors: - const auto xx1 = padded_vector_view(2, x.data().begin(), p_.lx); - const auto xx2 = padded_vector_view(2, x.data().begin() + p_.lx, p_.lx); - const auto xx3 = padded_vector_view(2, x.data().begin() + 2 * p_.lx, p_.lx); - const auto xx4 = padded_vector_view(2, x.data().begin() + 3 * p_.lx, p_.lx); - - // DO I = 1, N + 4 + 2 * kk - // laavec(I) = P(I) - // lacvec(I) = P(I + N + 4 + 2 * kk) - // muvec(I) = P(I + 2 * (N + 4 + 2 * kk)) - // gamvec(I) = P(I + 3 * (N + 4 + 2 * kk)) - // nn(I) = P(I + 4 * (N + 4 + 2 * kk)) - // ENDDO - // using views instead of vectors: + const auto xx1 = padded_vector_view<2>(x.data().begin(), p_.lx); + const auto xx2 = padded_vector_view<2>(x.data().begin() + p_.lx, p_.lx); + const auto xx3 = padded_vector_view<2>(x.data().begin() + 2 * p_.lx, p_.lx); + const auto xx4 = padded_vector_view<2>(x.data().begin() + 3 * p_.lx, p_.lx); + const auto chunk = p_.lx + 4 + 2 * p_.kk; const auto laavec = p_.P.data().begin(); const auto lacvec = p_.P.data().begin() + chunk; @@ -211,93 +130,54 @@ namespace { const auto gamvec = p_.P.data().begin() + 3 * chunk; const auto nn = p_.P.data().begin() + 4 * chunk; - // DO I = 3, N + 2 - // il1(I - 2) = I + kk - 1 - // il2(I - 2) = I + kk + 1 - // il3in3(I - 2) = I + kk - // il4(I - 2) = I + kk - 2 - // in1(I - 2) = I + 2 * kk - 1 - // in2ix2(I - 2) = I + 1 - // in4ix1(I - 2) = I - 1 - // ix3(I - 2) = I - // ix4(I - 2) = I - 2 - // ENDDO // using offsets into our views instead of vectors: - const int il1 = 2 + p_.kk - 1; - const int il2 = 2 + p_.kk + 1; - const int il3in3 = 2 + p_.kk; - const int il4 = 2 + p_.kk - 2; - const int in1 = 2 + 2 * p_.kk - 1; - const int in2ix2 = 2 + 1; // spilt in in2, ix2 at no cost? - const int in4ix1 = 2 - 1; // split in in4, ix1 at no cost? - const int ix3 = 2; - const int ix4 = 2 - 2; - - //DO I = 1, N - // FF1 = lacvec(il1(I)) * xx1(in4ix1(I)) - // FF1 = FF1 + laavec(il1(I) + 1) * xx2(in4ix1(I)) - // FF1 = FF1 + lacvec(il4(I) + 1) * xx2(ix4(I)) - // FF1 = FF1 + muvec(il2(I)) * nn(in2ix2(I)) * xx1(in2ix2(I)) - // FF1 = FF1 + muvec(il3in3(I) + 1) * xx2(ix3(I)) - // FFF = muvec(il3in3(I)) + lacvec(il3in3(I)) - // FF1 = FF1 - FFF * nn(il3in3(I)) * xx1(ix3(I)) - // dConc(I) = FF1 - gamvec(il3in3(I)) * xx1(ix3(I)) - // FF1 = gamvec(il3in3(I)) * xx1(ix3(I)) - // FF1 = FF1 + gamvec(il3in3(I)) * xx3(ix3(I)) - // FF1 = FF1 + gamvec(il3in3(I) + 1) * xx4(ix3(I)) - // FF1 = FF1 + lacvec(il1(I) + 1) * nn(in1(I)) * xx2(in4ix1(I)) - // FF1 = FF1 + muvec(il2(I) + 1) * nn(in2ix2(I)) * xx2(in2ix2(I)) - // FFF = muvec(il3in3(I) + 1) + lacvec(il3in3(I) + 1) - // FF1 = FF1 - FFF * nn(il3in3(I) + 1) * xx2(ix3(I)) - // FF1 = FF1 - laavec(il3in3(I) + 1) * xx2(ix3(I)) - // dConc(N + I) = FF1 - // FF1 = lacvec(il1(I)) * nn(in1(I)) * xx3(in4ix1(I)) - // FF1 = FF1 + laavec(il1(I) + 1) * xx4(in4ix1(I)) - // FF1 = FF1 + lacvec(il4(I) + 1) * xx4(ix4(I)) - // FF1 = FF1 + muvec(il2(I)) * nn(in2ix2(I)) * xx3(in2ix2(I)) - // FF1 = FF1 + muvec(il3in3(I) + 1) * xx4(ix3(I)) - // FFF = lacvec(il3in3(I)) + muvec(il3in3(I)) - // FF1 = FF1 - FFF * nn(il3in3(I)) * xx3(ix3(I)) - // FF1 = FF1 - gamvec(il3in3(I)) * xx3(ix3(I)) - // dConc(2 * N + I) = FF1 - // FF1 = lacvec(il1(I) + 1) * nn(in1(I)) * xx4(in4ix1(I)) - // FF1 = FF1 + muvec(il2(I) + 1) * nn(in2ix2(I)) * xx4(in2ix2(I)) - // FFF = lacvec(il3in3(I) + 1) + muvec(il3in3(I) + 1) - // FF1 = FF1 - FFF * nn(il3in3(I) + 1) * xx4(ix3(I)) - // FF1 = FF1 - gamvec(il3in3(I) + 1) * xx4(ix3(I)) - // dConc(3 * N + I) = FF1 - //ENDDO + const int nil2lx = 2; + const int il1 = nil2lx + p_.kk - 1; + const int il2 = nil2lx + p_.kk + 1; + const int il3 = nil2lx + p_.kk; + const int il4 = nil2lx + p_.kk - 2; + + const int in1 = nil2lx + 2 * p_.kk - 1; + const int in2 = nil2lx + 1; + const int in3 = nil2lx + p_.kk; + + const int ix1 = nil2lx - 1; + const int ix2 = nil2lx + 1; + const int ix3 = nil2lx; + const int ix4 = nil2lx - 2; + // using views into output vector: auto dx1 = dx.data().begin(); auto dx2 = dx1 + p_.lx; auto dx3 = dx2 + p_.lx; auto dx4 = dx3 + p_.lx; + for (int i = 0; i < p_.lx; ++i) { - dx1[i] = lacvec[il1 + i] * nn[in1 + i] * xx1[in4ix1 + i] - + laavec[il1 + i + 1] * xx2[in4ix1 + i] - + lacvec[il4 + i + 1] * xx2[ix4 + i] - + muvec[il2 + i] * nn[in2ix2 + i] * xx1[in2ix2 + i] - + muvec[il3in3 + i + 1] * xx2[ix3 + i] - - (muvec[il3in3 + i] + lacvec[il3in3 + i]) * nn[il3in3 + i] * xx1[ix3 + i] - - gamvec[il3in3 + i] * xx1[ix3 + i]; - dx2[i] = gamvec[il3in3 + i] * xx1[ix3 + i] - + gamvec[il3in3 + i] * xx3[ix3 + i] - + gamvec[il3in3 + i + 1] * xx4[ix3 + i] - + lacvec[il1 + i + 1] * nn[in1 + i] * xx2[in4ix1 + i] - + muvec[il2 + i + 1] * nn[in2ix2 + i] * xx2[in2ix2 + i] - - (muvec[il3in3 + 1 + i] + lacvec[il3in3 + 1 + i]) * nn[il3in3 + i + 1] * xx2[ix3 + i] - - laavec[il3in3 + i] * xx2[ix3 + i]; - dx3[i] = lacvec[il1 + i] * nn[in1 + i] * xx3[in4ix1 + i] - + laavec[il1 + i + 1] * xx4[in4ix1 + i] - + lacvec[il4 + i + 1] * xx4[ix4 + i] - + muvec[il2 + i] * nn[in2ix2 + i] * xx3[in2ix2 + i] - + muvec[il3in3 + i + 1] * xx4[ix3 + i] - - (lacvec[il3in3 + i] + muvec[il3in3 + i]) * nn[il3in3 + i] * xx3[ix3 + i] - - gamvec[il3in3 + i] * xx3[ix3 + i]; - dx4[i] = lacvec[il1 + i + 1] * nn[in1 + i] * xx4[in4ix1 + i] - + muvec[il2 + i + 1] * nn[in2ix2 + i] * xx4[in2ix2 + i] - - (lacvec[il3in3 + i + 1] + muvec[il3in3 + i + 1]) * nn[il3in3 + i + 1] * xx4[ix3 + i] - - gamvec[il3in3 + i + 1] * xx4[ix3 + i]; + dx1[i] = lacvec[il1 + i] * nn[in1 + i] * xx1[ix1 + i] + + laavec[il1 + i + 1] * xx2[ix1 + i] + + lacvec[il4 + i + 1] * xx2[ix4 + i] + + muvec[il2 + i] * nn[in2 + i] * xx1[ix2 + i] + + muvec[il3 + i + 1] * xx2[ix3 + i] + - (muvec[il3 + i] + lacvec[il3 + i]) * nn[in3 + i] * xx1[ix3 + i] + - gamvec[il3 + i] * xx1[ix3 + i]; + dx2[i] = gamvec[il3 + i] * xx1[ix3 + i] + + gamvec[il3 + i] * xx3[ix3 + i] + + gamvec[il3 + i + 1] * xx4[ix3 + i] + + lacvec[il1 + i + 1] * nn[in1 + i] * xx2[ix1 + i] + + muvec[il2 + i + 1] * nn[in2 + i] * xx2[ix2 + i] + - (muvec[il3 + 1 + i] + lacvec[il3 + 1 + i]) * nn[in3 + i + 1] * xx2[ix3 + i] + - laavec[il3 + i] * xx2[ix3 + i]; + dx3[i] = lacvec[il1 + i] * nn[in1 + i] * xx3[ix1 + i] + + laavec[il1 + i + 1] * xx4[ix1 + i] + + lacvec[il4 + i + 1] * xx4[ix4 + i] + + muvec[il2 + i] * nn[in2 + i] * xx3[ix2 + i] + + muvec[il3 + i + 1] * xx4[ix3 + i] + - (lacvec[il3 + i] + muvec[il3 + i]) * nn[in3 + i] * xx3[ix3 + i] + - gamvec[il3 + i] * xx3[ix3 + i]; + dx4[i] = lacvec[il1 + i + 1] * nn[in1 + i] * xx4[ix1 + i] + + muvec[il2 + i + 1] * nn[in2 + i] * xx4[ix2 + i] + - (lacvec[il3 + i + 1] + muvec[il3 + i + 1]) * nn[in3 + i + 1] * xx4[ix3 + i] + - gamvec[il3 + i + 1] * xx4[ix3 + i]; } } @@ -309,8 +189,10 @@ namespace { class cpp_daisie_cs_runmod_2 { public: - cpp_daisie_cs_runmod_2(param_t&& p) : - p_(p) + using jacobian = const_from_linear_rhs; + + cpp_daisie_cs_runmod_2(param_t&& par) : + p_(std::move(par)) { } @@ -319,22 +201,10 @@ namespace { { if (++p_.steps > max_cs_steps) throw std::runtime_error("cpp_daisie_cs_runmod_2: too many steps"); - // xx1 = c(0,0,x[1:lx],0) - // xx2 = c(0,0,x[(lx + 1):(2 * lx)],0) - // xx3 = c(0,0,x[(2 * lx + 1):(3 * lx)],0) - // using padded views instead of vectors: - const auto xx1 = padded_vector_view(2, x.data().begin(), p_.lx); - const auto xx2 = padded_vector_view(2, x.data().begin() + p_.lx, p_.lx); - const auto xx3 = padded_vector_view(2, x.data().begin() + 2 * p_.lx, p_.lx); - - // DO I = 1, N + 4 + 2 * kk - // laavec(I) = P(I) - // lacvec(I) = P(I + N + 4 + 2 * kk) - // muvec(I) = P(I + 2 * (N + 4 + 2 * kk)) - // gamvec(I) = P(I + 3 * (N + 4 + 2 * kk)) - // nn(I) = P(I + 4 * (N + 4 + 2 * kk)) - // ENDDO - // using views instead of vectors: + const auto xx1 = padded_vector_view<2>(x.data().begin(), p_.lx); + const auto xx2 = padded_vector_view<2>(x.data().begin() + p_.lx, p_.lx); + const auto xx3 = padded_vector_view<2>(x.data().begin() + 2 * p_.lx, p_.lx); + const auto chunk = p_.lx + 4 + 2 * p_.kk; const auto laavec = p_.P.data().begin(); const auto lacvec = p_.P.data().begin() + chunk; @@ -342,81 +212,47 @@ namespace { const auto gamvec = p_.P.data().begin() + 3 * chunk; const auto nn = p_.P.data().begin() + 4 * chunk; - // DO I = 3, N + 2 - // il1(I - 2) = I + kk - 1 - // il2(I - 2) = I + kk + 1 - // il3in3(I - 2) = I + kk - // il4(I - 2) = I + kk - 2 - // in1(I - 2) = I + 2 * kk - 1 - // in2ix2(I - 2) = I + 1 - // in4ix1(I - 2) = I - 1 - // ix3(I - 2) = I - // ix4(I - 2) = I - 2 - // ENDDO // using offsets into our views instead of vectors: - const int il1 = 2 + p_.kk - 1; - const int il2 = 2 + p_.kk + 1; - const int il3in3 = 2 + p_.kk; - const int il4 = 2 + p_.kk - 2; - const int in1 = 2 + 2 * p_.kk - 1; - const int in2ix2 = 2 + 1; // spilt in in2, ix2 at no cost? - const int in4ix1 = 2 - 1; // split in in4, ix1 at no cost? - const int ix3 = 2; - const int ix4 = 2 - 2; - - // DO I = 1, N - // FF1 = laavec(il1(I) + 1) * xx2(in4ix1(I)) - // FF1 = FF1 + lacvec(il4(I) + 1) * xx2(ix4(I)) - // FF1 = FF1 + muvec(il2(I) + 1) * xx2(ix3(I)) - // FF1 = FF1 + lacvec(il1(I)) * nn(in1(I)) * xx1(in4ix1(I)) - // FF1 = FF1 + muvec(il2(I)) * nn(in2ix2(I)) * xx1(in2ix2(I)) - // FFF = muvec(il3in3(I)) + lacvec(il3in3(I)) - // FF1 = FF1 - FFF * nn(il3in3(I)) * xx1(ix3(I)) - // FF1 = FF1 - gamvec(il3in3(I)) * xx1(ix3(I)) - // FFF = 0 - // IF(kk .EQ. 1) THEN - // FFF = laavec(il3in3(I)) * xx3(ix3(I)) - // FFF = FFF + 2 * lacvec(il1(I)) * xx3(in4ix1(I)) - // ENDIF - // dConc(I) = FF1 + FFF - // FF1 = gamvec(il3in3(I)) * xx1(ix3(I)) - // FF1 = FF1 + lacvec(il1(I) + 1) * nn(in1(I)) * xx2(in4ix1(I)) - // FF1 = FF1 + muvec(il2(I)+1) * nn(in2ix2(I)) * xx2(in2ix2(I)) - // FFF = muvec(il3in3(I) + 1) + lacvec(il3in3(I) + 1) - // FF1 = FF1 - FFF * nn(il3in3(I) + 1) * xx2(ix3(I)) - // FF1 = FF1 - laavec(il3in3(I) + 1) * xx2(ix3(I)) - // dConc(N + I) = FF1 - // FF1 = lacvec(il1(I)) * nn(in4ix1(I)) * xx3(in4ix1(I)) - // FF1 = FF1 + muvec(il2(I)) * nn(in2ix2(I)) * xx3(in2ix2(I)) - // FFF = lacvec(il3in3(I)) + muvec(il3in3(I)) - // FF1 = FF1 - FFF * nn(il3in3(I)) * xx3(ix3(I)) - // FF1 = FF1-(laavec(il3in3(I))+gamvec(il3in3(I)))*xx3(ix3(I)) - // dConc(2 * N + I) = FF1 - // ENDDO + const int nil2lx = 2; + const int il1 = nil2lx + p_.kk - 1; + const int il2 = nil2lx + p_.kk + 1; + const int il3 = nil2lx + p_.kk; + const int il4 = nil2lx + p_.kk - 2; + + const int in1 = nil2lx + 2 * p_.kk - 1; + const int in2 = nil2lx + 1; + const int in3 = nil2lx + p_.kk; + const int in4 = nil2lx-1; + + const int ix1 = nil2lx - 1; + const int ix2 = nil2lx + 1; + const int ix3 = nil2lx; + const int ix4 = nil2lx - 2; + // using views into output vector: auto dx1 = dx.data().begin(); auto dx2 = dx1 + p_.lx; auto dx3 = dx2 + p_.lx; + + const auto kk = (1 == p_.kk) ? 1.0 : 0.0; // make the loop body branch-free for (int i = 0; i < p_.lx; ++i) { - dx1[i] = laavec[il1 + i + 1] * xx2[in4ix1 + i] + dx1[i] = laavec[il1 + i + 1] * xx2[ix1 + i] + lacvec[il4 + i + 1] * xx2[ix4 + i] + muvec[il2 + i + 1] * xx2[ix3 + i] - + lacvec[il1 + i] * nn[in1 + i] * xx1[in4ix1 + i] - + muvec[il2 + i] * nn[in2ix2 + i] * xx1[in2ix2 + i] - - (muvec[il3in3 + i] + lacvec[il3in3 + i]) * nn[il3in3 + i] * xx1[ix3 + i] - - gamvec[il3in3 + i] * xx1[ix3 + i]; - if (1 == p_.kk) { - dx1[i] += laavec[il3in3 + i] * xx3[ix3 + i] + 2.0 * lacvec[il1 + i] * xx3[in4ix1 + i]; - } - dx2[i] = gamvec[il3in3 + i] * xx1[ix3 + i] - + lacvec[il1 + i + 1] * nn[in1 + i] * xx2[in4ix1 + i] - + muvec[il2 + i + 1] * nn[in2ix2 + i] * xx2[in2ix2 + i] - - (muvec[il3in3 + 1 + i] + lacvec[il3in3 + 1 + i]) * nn[il3in3 + i + 1] * xx2[ix3 + i] - - laavec[il3in3 + i] * xx2[ix3 + i]; - dx3[i] = lacvec[il1 + i] * nn[in4ix1 + i] * xx3[in4ix1 + i] - + muvec[il2 + i] * nn[in2ix2 + i] * xx3[in2ix2 + i] - - (lacvec[il3in3 + i] + muvec[il3in3 + i]) * nn[il3in3 + i] * xx3[ix3 + i] - - (laavec[il3in3 + i] + gamvec[il3in3 + i]) * xx3[ix3 + i]; + + lacvec[il1 + i] * nn[in1 + i] * xx1[ix1 + i] + + muvec[il2 + i] * nn[in2 + i] * xx1[ix2 + i] + - (muvec[il3 + i] + lacvec[il3 + i]) * nn[in3 + i] * xx1[ix3 + i] + - gamvec[il3 + i] * xx1[ix3 + i] + + kk * (laavec[il3 + i] * xx3[ix3 + i] + 2.0 * lacvec[il1 + i] * xx3[ix1 + i]); + dx2[i] = gamvec[il3 + i] * xx1[ix3 + i] + + lacvec[il1 + i + 1] * nn[in1 + i] * xx2[ix1 + i] + + muvec[il2 + i + 1] * nn[in2 + i] * xx2[ix2 + i] + - (muvec[il3 + 1 + i] + lacvec[il3 + 1 + i]) * nn[in3 + i + 1] * xx2[ix3 + i] + - laavec[il3 + i] * xx2[ix3 + i]; + dx3[i] = lacvec[il1 + i] * nn[in4 + i] * xx3[ix1 + i] + + muvec[il2 + i] * nn[in2 + i] * xx3[ix2 + i] + - (lacvec[il3 + i] + muvec[il3 + i]) * nn[in3 + i] * xx3[ix3 + i] + - (laavec[il3 + i] + gamvec[il3 + i]) * xx3[ix3 + i]; } } @@ -430,7 +266,7 @@ namespace { //' Driver for the boost::odeint solver //' //' @name daisie_odeint_cs -RcppExport SEXP daisie_odeint_cs(SEXP rrunmod, SEXP ry, SEXP rtimes, SEXP rlx, SEXP rkk, SEXP rpar, SEXP Stepper, SEXP atolint, SEXP reltolint) { +RcppExport SEXP daisie_odeint_cs(SEXP rrunmod, SEXP ry, SEXP rtimes, SEXP rlx, SEXP rkk, SEXP rpar, SEXP Stepper, SEXP ratol, SEXP rrtol) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; @@ -440,8 +276,8 @@ BEGIN_RCPP auto lx = as(rlx); auto kk = as(rkk); auto stepper = as(Stepper); - auto atol = as(atolint); - auto rtol = as(reltolint); + auto atol = as(ratol); + auto rtol = as(rrtol); auto p = param_t(lx, kk, as(rpar)); if (runmod == "daisie_runmod") { diff --git a/src/DAISIE_IW.cpp b/src/DAISIE_IW.cpp index e6f55b53..27d39945 100644 --- a/src/DAISIE_IW.cpp +++ b/src/DAISIE_IW.cpp @@ -19,6 +19,7 @@ using namespace Eigen; // num_threads unsigned daisie_odeint_iw_num_threads_ = std::max(1u, std::thread::hardware_concurrency()); +using namespace daisie_odeint::jacobian_policy; namespace { @@ -178,6 +179,8 @@ namespace { struct daisie_iw_wrapper { + using jacobian = const_from_linear_rhs; + std::unique_ptr pool; std::unique_ptr dev; diff --git a/src/DAISIE_odeint.h b/src/DAISIE_odeint.h index 4bea02d3..af595928 100644 --- a/src/DAISIE_odeint.h +++ b/src/DAISIE_odeint.h @@ -1,124 +1,113 @@ +#pragma once #ifndef DAISIE_ODEINT_H_INCLUDED #define DAISIE_ODEINT_H_INCLUDED -// boiler-plate code calling into boost::odeint - - -#define STRICT_R_HEADERS -#include - -// [[Rcpp::plugins(cpp14)]] -// [[Rcpp::depends(BH)]] - -#include - -// Provide Forward Declarations -namespace Rcpp { - - namespace traits{ - - // Setup non-intrusive extension via template specialization for - // 'ublas' class boost::numeric::ublas - - // Support for wrap - template SEXP wrap(const boost::numeric::ublas::vector & obj); - - // Support for as - template class Exporter< boost::numeric::ublas::vector >; - - } -} - - -#include +#include "ublas_types.h" #include -#include #include #include -// boost::numeric::ublas wrapping from: -// https://gallery.rcpp.org/articles/custom-templated-wrap-and-as-for-seamingless-interfaces/ -namespace Rcpp { - - namespace traits{ - - // Defined wrap case - template SEXP wrap(const boost::numeric::ublas::vector & obj){ - const int RTYPE = Rcpp::traits::r_sexptype_traits::rtype ; - - return Rcpp::Vector< RTYPE >(obj.begin(), obj.end()); - }; +using namespace Rcpp; +using namespace boost::numeric::odeint; - // Defined as< > case - template class Exporter< boost::numeric::ublas::vector > { - typedef typename boost::numeric::ublas::vector OUT ; +// type of the ode state +using state_type = vector_t; - // Convert the type to a valid rtype. - const static int RTYPE = Rcpp::traits::r_sexptype_traits< T >::rtype ; - Rcpp::Vector vec; - public: - Exporter(SEXP x) : vec(x) { - if (TYPEOF(x) != RTYPE) - throw std::invalid_argument("Wrong R type for mapped 1D array"); - } - OUT get() { - // Need to figure out a way to perhaps do a pointer pass? - OUT x(vec.size()); - std::copy(vec.begin(), vec.end(), x.begin()); // have to copy data - return x; - } - }; +// zero-value padded view into vector +template +class padded_vector_view +{ +public: + padded_vector_view(const double* data, int n) : + sdata_(data - Pad), sn_(n + Pad) + { + } + // returns 0.0 for indices 'i' outside [Pad, Pad + n) + double operator[](int i) const + { + return (i >= Pad && i < sn_) ? *(sdata_ + i) : 0.0; } -} +private: + const double* sdata_ = nullptr; // sdata_[Pad] == data[0] + const int sn_ = 0; +}; -using namespace Rcpp; -using namespace boost::numeric::odeint; +namespace daisie_odeint { -// type of the ode state -using state_type = boost::numeric::ublas::vector; + extern double abm_factor; + + template + inline void do_integrate(double atol, double rtol, Rhs rhs, state_type& y, double t0, double t1) + { + integrate_adaptive(make_controlled(atol, rtol), rhs, y, t0, t1, 0.1 * (t1 - t0)); + } -namespace daisie_odeint { + template + inline void abm(Rhs rhs, state_type& y, double t0, double t1) + { + auto abm = adams_bashforth_moulton{}; + abm.initialize(rhs, y, t0, abm_factor * (t1 - t0)); + integrate_const(abm, rhs, y, t0, t1, abm_factor * (t1 - t0)); + } -extern double abm_factor; + template + inline void ab(Rhs rhs, state_type& y, double t0, double t1) + { + auto ab = adams_bashforth{}; + ab.initialize(rhs, y, t0, abm_factor * (t1 - t0)); + integrate_const(ab, rhs, y, t0, t1, abm_factor * (t1 - t0)); + } -template -inline void do_integrate(double atol, double rtol, Rhs rhs, state_type& y, double t0, double t1) -{ - integrate_adaptive(make_controlled(atol, rtol), rhs, y, t0, t1, 0.1 * (t1 - t0)); -} + namespace jacobian_policy { + // Evaluator of the Jacobian for linear, time independent systems + // dxdt = Ax => Jacobian = t(A) + template + struct const_from_linear_rhs + { + explicit const_from_linear_rhs(RHS& rhs) : rhs_(rhs) + { + } -template -inline void abm(Rhs rhs, state_type& y, double t0, double t1) -{ - auto abm = adams_bashforth_moulton{}; - abm.initialize(rhs, y, t0, abm_factor * (t1 - t0)); - integrate_const(abm, rhs, y, t0, t1, abm_factor * (t1 - t0)); -} + void operator()(const vector_t& x, matrix_t& J, double t, vector_t& /*dfdt*/) + { + if (!J_) { + // once-only, generic evaluation + J_ = std::make_unique>(J.size1(), J.size2()); + auto single = vector_t(x.size(), 0); + auto dxdt = vector_t(x.size()); + for (size_t i = 0; i < J.size1(); ++i) { + single[i] = 1.0; + auto col = ublas::matrix_column>(*J_, i); + std::copy(col.begin(), col.end(), dxdt.begin()); + rhs_(single, dxdt, 0); + std::copy(dxdt.begin(), dxdt.end(), col.begin()); + single[i] = 0.0; + } + } + J = *J_; + } + RHS& rhs_; + std::unique_ptr> J_; + }; -template -inline void ab(Rhs rhs, state_type& y, double t0, double t1) -{ - auto ab = adams_bashforth{}; - ab.initialize(rhs, y, t0, abm_factor * (t1 - t0)); - integrate_const(ab, rhs, y, t0, t1, abm_factor * (t1 - t0)); -} + } -// wrapper around odeint::integrate + // wrapper around odeint::integrate // maps runtime stepper name -> compiletime odeint::stepper type template inline void integrate( @@ -155,7 +144,7 @@ inline void ab(Rhs rhs, state_type& y, double t0, double t1) case '6': ab<6>(rhs, y, t0, t1); break; case '7': ab<7>(rhs, y, t0, t1); break; case '8': ab<8>(rhs, y, t0, t1); break; - default: throw std::runtime_error("DAISIE_odeint_helper::integrate: unsupported steps for admam_bashforth_moulton"); + default: throw std::runtime_error("DAISIE_odeint_helper::integrate: unsupported steps for admam_bashforth"); } } else if (0 == stepper.compare(0, stepper.size() - 2, "odeint::adams_bashforth_moulton")) { @@ -172,6 +161,14 @@ inline void ab(Rhs rhs, state_type& y, double t0, double t1) default: throw std::runtime_error("DAISIE_odeint_helper::integrate: unsupported steps for admam_bashforth_moulton"); } } + else if ("odeint::rosenbrock4" == stepper) { + // another outlier in calling convention + using stepper_t = rosenbrock4; + using controlled_stepper_t = rosenbrock4_controller; + auto jac = typename Rhs::type::jacobian(rhs); + auto sys = std::make_pair(std::ref(rhs), std::ref(jac)); + integrate_adaptive(controlled_stepper_t(atol, rtol), sys, y, t0, t1, 0.1 * (t1 - t0)); + } else { throw std::runtime_error("DAISIE_odeint_helper::integrate: unknown stepper"); } diff --git a/src/DAISIE_types.h b/src/DAISIE_types.h new file mode 100644 index 00000000..1f21d261 --- /dev/null +++ b/src/DAISIE_types.h @@ -0,0 +1,2 @@ +#pragma once +#include "odeint_types.h" diff --git a/src/ublas_types.h b/src/ublas_types.h new file mode 100644 index 00000000..125cc936 --- /dev/null +++ b/src/ublas_types.h @@ -0,0 +1,116 @@ +#pragma once +#ifndef UBLAS_TYPES_H_INCLUDED +#define UBLAS_TYPES_H_INCLUDED + +#include +#include +#include + + +namespace ublas = boost::numeric::ublas; + + +template using vector_t = ublas::vector; +template using matrix_t = ublas::matrix; + + +// forward declarations Rcpp <-> boost::numeric::ublas +namespace Rcpp { + + namespace traits { + + template SEXP wrap(const vector_t&); + template SEXP wrap(const matrix_t&); + + template vector_t as(SEXP); + template matrix_t as(SEXP); + + template class Exporter>; + template class Exporter>; + + } + +} + + +#include + + +namespace Rcpp { + + namespace traits { + + template inline SEXP wrap(const vector_t& obj) { + const int RTYPE = Rcpp::traits::r_sexptype_traits::rtype; + return Rcpp::Vector(obj.begin(), obj.end()); + } + + + template inline SEXP wrap(const matrix_t& obj) { + const size_t nr = static_cast(obj.size1()); + const size_t nc = static_cast(obj.size2()); + const int RTYPE = Rcpp::traits::r_sexptype_traits::rtype; + Rcpp::Matrix rmat(nr, nc); + for (size_t i = 0; i < nr; ++i) { + for (size_t j = 0; j < nc; ++j) { + rmat(i, j) = obj(i, j); + } + } + return rmat; + } + + + template + class Exporter> + { + private: + static constexpr int RTYPE = Rcpp::traits::r_sexptype_traits::rtype; + Rcpp::Vector rvec; + + public: + Exporter(SEXP x) : rvec(x) { + if (TYPEOF(x) != RTYPE) { + throw std::invalid_argument("Wrong R type for mapped 1D array"); + } + } + + vector_t get() { + vector_t x(rvec.size()); + std::copy(rvec.begin(), rvec.end(), x.begin()); + return x; + } + }; + + + template + class Exporter> + { + private: + static constexpr int RTYPE = Rcpp::traits::r_sexptype_traits::rtype; + Rcpp::Matrix rmat; + + public: + Exporter(SEXP x) : rmat(x) { + if (TYPEOF(x) != RTYPE) { + throw std::invalid_argument("Wrong R type for mapped 2D array"); + } + } + + matrix_t get() { + const size_t nr = static_cast(rmat.rows()); + const size_t nc = static_cast(rmat.cols()); + matrix_t x(nr, nc); + for (size_t i = 0; i < nr; ++i) { + for (size_t j = 0; j < nc; ++j) { + x(i, j) = rmat(i, j); + } + } + return x; + } + }; + + } + +} + +#endif diff --git a/tests/testthat/test-DAISIE_ML4.R b/tests/testthat/test-DAISIE_ML4.R index b6642552..9449ba2a 100644 --- a/tests/testthat/test-DAISIE_ML4.R +++ b/tests/testthat/test-DAISIE_ML4.R @@ -3,12 +3,16 @@ test_that("DAISIE_ML4 is silent and produces correct output", { utils::data(Galapagos_datalist) loglik <- DAISIE_ML4( datalist = Galapagos_datalist, - initparsopt = c(2.5, 2.7, 20, 0.009, 1.01, 2), + initparsopt = c(1.05, 0.36, 26.6, 0.0029, 0.73, 0.1), idparsopt = 1:6, parsfix = NULL, idparsfix = NULL, + methode = 'lsodes', + optimmethod = 'simplex', CS_version = create_CS_version(model = 2, - relaxed_par = "cladogenesis")) + relaxed_par = "extinction", + par_sd = 0.1, + par_upper_bound = 1)) expect_equal(2, 2) }) @@ -27,7 +31,9 @@ test_that("DAISIE_loglik_all_choosepar4 is silent and produces correct output", datalist = Galapagos_datalist, methode = "lsodes", CS_version = list(model = 2, - relaxed_par = "cladogenesis"), + relaxed_par = "cladogenesis", + par_sd = 1, + par_upper_bound = Inf), abstolint = 1e-16, reltolint = 1e-10 )))) diff --git a/tests/testthat/test-DAISIE_format_IW.R b/tests/testthat/test-DAISIE_format_IW.R index b2506e96..246d0622 100644 --- a/tests/testthat/test-DAISIE_format_IW.R +++ b/tests/testthat/test-DAISIE_format_IW.R @@ -540,24 +540,24 @@ test_that("silent when species with two trait states with brts_table <- matrix(ncol = 5, nrow = 4) colnames(brts_table) <- c("brt", "clade", "event", "endemic", "col") brts_table[1, ] <- c(5.00000000000000, 0, 0, NA, NA) - brts_table[2, ] <- c(3.10634202528338, 1, 1, 1, NA) - brts_table[3, ] <- c(1.52330128016821, 2, 1, 1, NA) - brts_table[4, ] <- c(1.28012784155125, 2, 2, 1, NA) + brts_table[2, ] <- c(3.10261367452990, 1, 1, 1, NA) + brts_table[3, ] <- c(1.50562999775257, 2, 1, 1, NA) + brts_table[4, ] <- c(1.26245655913561, 2, 2, 1, NA) expected_IW_format[[1]][[1]] <- list(island_age = 5, not_present = 13, stt_all = stt_all, brts_table = brts_table) expected_IW_format[[1]][[2]] <- list( - branching_times = c(5.00000000000000, - 3.10634202528338), + branching_times = c(5.0000000000000, + 3.1026136745299), stac = 2, missing_species = 0 ) expected_IW_format[[1]][[3]] <- list( branching_times = c(5.00000000000000, - 1.52330128016821, - 1.28012784155125), + 1.50562999775257, + 1.26245655913561), stac = 2, missing_species = 0 ) diff --git a/tests/testthat/test-DAISIE_loglik_CS.R b/tests/testthat/test-DAISIE_loglik_CS.R index 071c5036..1de5544d 100644 --- a/tests/testthat/test-DAISIE_loglik_CS.R +++ b/tests/testthat/test-DAISIE_loglik_CS.R @@ -30,7 +30,8 @@ test_that("DAISIE_loglik_CS_choice produces correct output for relaxed-rate missnumspec <- 0 CS_version <- list(model = 2, relaxed_par = "cladogenesis", - sd = 1) + par_sd = 1, + par_upper_bound = Inf) invisible(capture.output(loglik <- DAISIE_loglik_CS_choice(pars1 = pars1, pars2 = pars2, @@ -85,7 +86,8 @@ test_that("DAISIE_loglik_all produces correct output for relaxed-rate model", { methode = "lsodes", CS_version = list(model = 2, relaxed_par = "cladogenesis", - sd = 1), + par_sd = 1, + par_upper_bound = Inf), abstolint = 1e-16, reltolint = 1e-10 ) diff --git a/tests/testthat/test-DAISIE_loglik_integrate.R b/tests/testthat/test-DAISIE_loglik_integrate.R index 7b7cdf18..1aad41c9 100644 --- a/tests/testthat/test-DAISIE_loglik_integrate.R +++ b/tests/testthat/test-DAISIE_loglik_integrate.R @@ -8,7 +8,8 @@ test_that("DAISIE_loglik_integrate produces correct ouput on single lineage", { methode <- "lsodes" CS_version <- list(model = 2, relaxed_par = 'carrying_capacity', - sd = 2) + par_sd = 2, + par_upper_bound = Inf) abstolint <- 1e-16 reltolint <- 1e-10 verbose <- FALSE @@ -38,7 +39,8 @@ test_that("DAISIE_loglik_integrate produces correct ouput on radiation", { methode <- "lsodes" CS_version <- list(model = 2, relaxed_par = 'carrying_capacity', - sd = 10) + par_sd = 10, + par_upper_bound = Inf) abstolint <- 1e-16 reltolint <- 1e-10 verbose <- FALSE diff --git a/tests/testthat/test-DAISIE_utils.R b/tests/testthat/test-DAISIE_utils.R index f954bdb7..f3d90d6d 100644 --- a/tests/testthat/test-DAISIE_utils.R +++ b/tests/testthat/test-DAISIE_utils.R @@ -232,15 +232,23 @@ test_that("create_CS_version produces correct output", { CS_version <- create_CS_version(model = 1, relaxed_par = NULL) expect_equal(CS_version, list(model = 1, - relaxed_par = NULL)) + relaxed_par = NULL, + par_sd = 0, + par_upper_bound = Inf)) CS_version <- create_CS_version(model = 2, - relaxed_par = "cladogenesis") + relaxed_par = "cladogenesis", + par_sd = 10, + par_upper_bound = Inf) expect_equal(CS_version, list(model = 2, - relaxed_par = "cladogenesis")) + relaxed_par = "cladogenesis", + par_sd = 10, + par_upper_bound = Inf)) CS_version <- create_CS_version(model = 3, relaxed_par = NULL) expect_equal(CS_version, list(model = 3, - relaxed_par = NULL)) + relaxed_par = NULL, + par_sd = 0, + par_upper_bound = Inf)) }) diff --git a/tests/testthat/test_odeint.R b/tests/testthat/test_odeint.R new file mode 100644 index 00000000..6a47f4a9 --- /dev/null +++ b/tests/testthat/test_odeint.R @@ -0,0 +1,65 @@ +test_that("odeint solvers give the same result as deSolve solvers", { + utils::data(Galapagos_datalist_2types) + pars1 <- c( + 0.195442017, + 0.087959583, + Inf, + 0.002247364, + 0.873605049, + 3755.202241, + 8.909285094, + 14.99999923, + 0.002247364, + 0.873605049, + 0.163 + ) + pars2 <- c(40, 11, 0, 1) + methode <- 'deSolve_R::lsodes' + loglik_lsodes_R <- DAISIE_loglik_all(pars1, pars2, Galapagos_datalist_2types, methode = methode) + methode <- 'deSolve_R::lsoda' + loglik_lsoda_R <- DAISIE_loglik_all(pars1, pars2, Galapagos_datalist_2types, methode = methode) + methode <- 'lsodes' + loglik_lsodes <- DAISIE_loglik_all(pars1, pars2, Galapagos_datalist_2types, methode = methode) + methode <- 'lsoda' + loglik_lsoda <- DAISIE_loglik_all(pars1, pars2, Galapagos_datalist_2types, methode = methode) + expect_equal(loglik_lsoda_R,loglik_lsoda) + expect_equal(loglik_lsodes_R,loglik_lsodes) + expect_equal(loglik_lsoda,loglik_lsodes) + methode <- 'odeint::runge_kutta_cash_karp54' + loglik_rkck54 <- DAISIE_loglik_all(pars1, pars2, Galapagos_datalist_2types, methode = methode) + methode <- 'odeint::runge_kutta_fehlberg78' + loglik_rkf78 <- DAISIE_loglik_all(pars1, pars2, Galapagos_datalist_2types, methode = methode) + methode <- 'odeint::runge_kutta_dopri5' + loglik_rkd5 <- DAISIE_loglik_all(pars1, pars2, Galapagos_datalist_2types, methode = methode) + methode <- 'odeint::bulirsch_stoer' + loglik_bs <- DAISIE_loglik_all(pars1, pars2, Galapagos_datalist_2types, methode = methode) + methode <- 'odeint::rosenbrock4' + loglik_rb <- DAISIE_loglik_all(pars1, pars2, Galapagos_datalist_2types, methode = methode) + methode <- 'odeint::adams_bashforth_moulton_1' + DAISIE_CS_max_steps(100000000) + DAISIE_abm_factor(0.000001) + loglik_abm <- DAISIE_loglik_all(pars1, pars2, Galapagos_datalist_2types, methode = methode) + expect_equal(loglik_lsodes, loglik_rkck54) + expect_equal(loglik_lsodes, loglik_rkf78) + expect_equal(loglik_lsodes, loglik_rkd5) + expect_equal(loglik_lsodes, loglik_bs) + expect_equal(loglik_lsodes, loglik_rb) + expect_equal(loglik_lsodes, loglik_abm, tol = 1E-6) + + pars1a <- pars1 + pars1a[6] <- Inf + methode <- 'deSolve_R::lsoda' + loglik_lsoda_R_Inf <- DAISIE_loglik_all(pars1a, pars2, Galapagos_datalist_2types, methode = methode) + methode <- 'lsoda' + loglik_lsoda_F_Inf <- DAISIE_loglik_all(pars1a, pars2, Galapagos_datalist_2types, methode = methode) + expect_equal(loglik_lsoda_R_Inf,loglik_lsoda_F_Inf) + + methode <- 'deSolve_R::lsodes' + loglik_lsodes_R_Inf <- DAISIE_loglik_all(pars1a, pars2, Galapagos_datalist_2types, methode = methode) + methode <- 'lsodes' + loglik_lsodes_F_Inf <- DAISIE_loglik_all(pars1a, pars2, Galapagos_datalist_2types, methode = methode) + expect_equal(loglik_lsodes_R_Inf,loglik_lsodes_F_Inf) + + expect_equal(loglik_lsoda_R_Inf,loglik_lsoda_R, tol = 1E-4) + expect_equal(loglik_lsodes_R_Inf,loglik_lsodes_R, tol = 1E-4) +}) diff --git a/tests/testthat/test_printing.R b/tests/testthat/test_printing.R new file mode 100644 index 00000000..4d331126 --- /dev/null +++ b/tests/testthat/test_printing.R @@ -0,0 +1,8 @@ +test_that("printing is done", { + print_parameters_and_loglik(pars = c(1,2:6), loglik = -3, verbose = T, type = 'clade_loglik') + print_parameters_and_loglik(pars = c(2:6), loglik = -3, verbose = T, type = 'island_loglik') + print_parameters_and_loglik(pars = c(2:6), loglik = -3, verbose = T, type = 'island_ML') + print_parameters_and_loglik(pars = c(1:11), loglik = -3, verbose = T, type = 'global_ML', distance_dep = 'power') + print_parameters_and_loglik(pars = c(1:11), loglik = -3, verbose = T, type = 'global_ML', distance_dep = 'sigmoidal_col') + print_parameters_and_loglik(pars = data.frame(rbind(c(2:6),c(12:16))), loglik = -3, verbose = T, type = 'multiple_island_ML') +})