diff --git a/.github/workflows/R-CMD-check.yaml b/.github/workflows/R-CMD-check.yaml index e96e2b2c..569d8a9d 100644 --- a/.github/workflows/R-CMD-check.yaml +++ b/.github/workflows/R-CMD-check.yaml @@ -6,6 +6,8 @@ on: - main - master - develop + - DAISIE_IW_B + pull_request: branches: - main @@ -84,13 +86,6 @@ jobs: run: find check -name 'testthat.Rout*' -exec cat '{}' \; || true shell: bash - - name: Install covr and test coverage - if: matrix.config.r == 'devel' - run: | - remotes::install_cran("covr") - covr::codecov() - shell: Rscript {0} - - name: Upload check results if: failure() uses: actions/upload-artifact@main diff --git a/.github/workflows/test-coverage.yaml b/.github/workflows/test-coverage.yaml new file mode 100644 index 00000000..2ca07e37 --- /dev/null +++ b/.github/workflows/test-coverage.yaml @@ -0,0 +1,51 @@ +on: + push: + branches: + - main + - master + - develop + pull_request: + branches: + - main + - master + - develop + +name: test-coverage + +jobs: + test-coverage: + runs-on: ubuntu-20.04 + + env: + GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} + steps: + - uses: actions/checkout@v2 + + - uses: r-lib/actions/setup-r@v1 + + - uses: r-lib/actions/setup-pandoc@v1 + + - name: Query dependencies + run: | + install.packages('remotes') + saveRDS(remotes::dev_package_deps(dependencies = TRUE), ".github/depends.Rds", version = 2) + writeLines(sprintf("R-%i.%i", getRversion()$major, getRversion()$minor), ".github/R-version") + shell: Rscript {0} + + - name: Restore R package cache + uses: actions/cache@v2 + with: + path: ${{ env.R_LIBS_USER }} + key: ${{ runner.os }}-${{ hashFiles('.github/R-version') }}-1-${{ hashFiles('.github/depends.Rds') }} + restore-keys: ${{ runner.os }}-${{ hashFiles('.github/R-version') }}-1- + + - name: Install dependencies + run: | + install.packages(c("remotes")) + remotes::install_deps(dependencies = TRUE) + remotes::install_cran("covr") + shell: Rscript {0} + + - name: Test coverage + run: covr::codecov() + shell: Rscript {0} diff --git a/DESCRIPTION b/DESCRIPTION index 5d4a79b1..1bfbd45d 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,8 +1,8 @@ Package: DAISIE Type: Package Title: Dynamical Assembly of Islands by Speciation, Immigration and Extinction -Version: 3.2.1 -Date: 2021-03-23 +Version: 4.0.2 +Date: 2021-05-26 Depends: R (>= 3.5.0) biocViews: SystemRequirements: C++14 @@ -34,7 +34,7 @@ Suggests: gridExtra, dplyr, ggplot2, - ggtree, + ggtree (>= 3.0.0), tidytree, tidyr, purrr, diff --git a/NEWS.md b/NEWS.md index aa05dbe5..444942fa 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,14 +1,38 @@ +# DAISIE 4.0.2 + +**N.B.: MacOS users may experience issues when installing DAISIE, especially when on MacOS Big Sur. If that is you case, please see [here](https://github.com/rsetienne/DAISIE/blob/6da0e3f65680d5f237345ef80935bda7541cf230/doc/DAISIE_macOS.md) for detailed installation instructions.** + +* Suggest ggtree >= 3.0.0. + +# DAISIE 4.0.1 + +**N.B.: MacOS users may experience issues when installing DAISIE, especially when on MacOS Big Sur. If that is you case, please see [here](https://github.com/rsetienne/DAISIE/blob/6da0e3f65680d5f237345ef80935bda7541cf230/doc/DAISIE_macOS.md) for detailed installation instructions.** + +* Fix possibility of fitting CS model with IW likelihood on simulated data by setting `CS_version = 0`. Improve `CS_version` documentation. + +# DAISIE 4.0.0 + +**N.B.: MacOS users may experience issues when installing DAISIE, especially when on MacOS Big Sur. If that is you case, please see [here](https://github.com/rsetienne/DAISIE/blob/6da0e3f65680d5f237345ef80935bda7541cf230/doc/DAISIE_macOS.md) for detailed installation instructions.** + +* Fix bug when calculating conditional probabilities, which are now correctly calculated from island age to the present. +* Fix bug when calculating probabilities upon migration, which assumed no recolonisation was possible in the CS model. Handling recolonisation in the same manner is not possible for the IW model, so an approximation is now made. The influence of the bug and of the approximation in the IW model is expected to be minimal, particularly in cases where the colonisation rate is low. Such cases of low colonisation are the norm. +* Add `num_cycles` argument to ML functions, allowing to specify how many cycles the optimizer should take. Defaults to 1. +* Minor vignette corrections. +* Improve README.md documentation. + # DAISIE 3.2.1 +[![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.4633441.svg)](https://doi.org/10.5281/zenodo.4633441) -**N.B.: MacOS users may experience issues when installing DAISIE, especially when on MacOS Big Sur. If that is you case, please see [here](doc/DAISIE_macOS.md) for detailed installation instructions.** +**N.B.: MacOS users may experience issues when installing DAISIE, especially when on MacOS Big Sur. If that is you case, please see [here](https://github.com/rsetienne/DAISIE/blob/6da0e3f65680d5f237345ef80935bda7541cf230/doc/DAISIE_macOS.md) for detailed installation instructions.** * Minor documentation improvements. # DAISIE 3.2.0 +[![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.4630604.svg)](https://doi.org/10.5281/zenodo.4630604) + -**N.B.: MacOS users may experience issues when installing DAISIE, especially when on MacOS Big Sur. If that is you case, please see [here](doc/DAISIE_macOS.md) for detailed installation instructions.** +**N.B.: MacOS users may experience issues when installing DAISIE, especially when on MacOS Big Sur. If that is you case, please see [here](https://github.com/rsetienne/DAISIE/blob/6da0e3f65680d5f237345ef80935bda7541cf230/doc/DAISIE_macOS.md) for detailed installation instructions.** -## Changes * `DAISIE_loglikg_IW()` is now more efficient and numerically stable. Numerical integration is now done via C++ with package `odeint`. * Add relaxed rate capabilities (both inference and simulations). Relaxed rate models allow for parameters to not be static, but to be sampled by specific probability distributions. * Introduce `MinAge` data status in DAISIE data objects. A status containing `MinAge` sets a lower boundary for colonization in situations when the precise colonization time is unknown. This is interpreted by `DAISIE_dataprep()` so that the information is passed on to the likelihood optimization functions. See the `DAISIE_dataprep()` help page for more details. In the back-end this results in new `stac` values 8 and 9. diff --git a/R/DAISIE_ExpEIN.R b/R/DAISIE_ExpEIN.R index e5a1ca59..aed3d5ea 100644 --- a/R/DAISIE_ExpEIN.R +++ b/R/DAISIE_ExpEIN.R @@ -45,7 +45,7 @@ DAISIE_ExpEIN <- function(t, pars, M, initEI = c(0, 0)) { DD <- laa + 2 * lac E0 <- initEI[1] I0 <- initEI[2] - if (t == Inf) { + if (t[1] == Inf) { Imm <- ga * M2 / B End <- DD / A * Imm } else { diff --git a/R/DAISIE_ML1.R b/R/DAISIE_ML1.R index fd76f70c..6b0dad08 100644 --- a/R/DAISIE_ML1.R +++ b/R/DAISIE_ML1.R @@ -95,7 +95,8 @@ DAISIE_ML1 <- function( verbose = 0, tolint = c(1E-16, 1E-10), island_ontogeny = NA, - jitter = 0) { + jitter = 0, + num_cycles = 1) { # datalist = list of all data: branching times, status of clade, and numnber of missing species # datalist[[,]][1] = list of branching times (positive, from present to past) # - max(brts) = age of the island @@ -297,7 +298,8 @@ DAISIE_ML1 <- function( CS_version = CS_version, abstolint = tolint[1], reltolint = tolint[2], - jitter = jitter + jitter = jitter, + num_cycles = num_cycles ) if (out$conv != 0) { cat( diff --git a/R/DAISIE_ML2.R b/R/DAISIE_ML2.R index 0fc72785..3cff782b 100644 --- a/R/DAISIE_ML2.R +++ b/R/DAISIE_ML2.R @@ -77,7 +77,8 @@ DAISIE_ML2 <- function( optimmethod = "subplex", verbose = 0, tolint = c(1E-16, 1E-10), - jitter = 0) { + jitter = 0, + num_cycles = 1) { # datalist = list of all data: branching times, status of clade, and numnber of missing species # datalist[[,]][1] = list of branching times (positive, from present to past) # - max(brts) = age of the island @@ -169,7 +170,8 @@ DAISIE_ML2 <- function( methode = methode, abstolint = tolint[1], reltolint = tolint[2], - jitter = jitter + jitter = jitter, + num_cycles = num_cycles ) if (out$conv != 0) { cat("Optimization has not converged. Try again with different initial values.\n") diff --git a/R/DAISIE_ML3.R b/R/DAISIE_ML3.R index 1e71cd70..0017b46c 100644 --- a/R/DAISIE_ML3.R +++ b/R/DAISIE_ML3.R @@ -67,7 +67,8 @@ DAISIE_ML3 <- function( CS_version = 1, verbose = 0, tolint = c(1E-16, 1E-10), - jitter = 0) { + jitter = 0, + num_cycles = 1) { # datalist = list of all data: branching times, status of clade, and numnber of missing species # datalist[[,]][1] = list of branching times (positive, from present to past) # - max(brts) = age of the island @@ -166,7 +167,8 @@ DAISIE_ML3 <- function( CS_version = CS_version, abstolint = tolint[1], reltolint = tolint[2], - jitter = jitter + jitter = jitter, + num_cycles = num_cycles ) if (out$conv != 0) { cat("Optimization has not converged. Try again with different initial values.\n") diff --git a/R/DAISIE_ML4.R b/R/DAISIE_ML4.R index a14af291..a08cbbe1 100644 --- a/R/DAISIE_ML4.R +++ b/R/DAISIE_ML4.R @@ -64,7 +64,7 @@ DAISIE_loglik_all_choosepar4 <- function(trparsopt, #' \item{loglik}{ gives the maximum loglikelihood} #' \item{df}{ #' gives the number -#' of estimated parameters, i.e. degrees of feedom} +#' of estimated parameters, i.e. degrees of freedom} #' \item{conv}{ gives a #' message on convergence of optimization; conv = 0 means convergence} #' @keywords internal @@ -86,7 +86,8 @@ DAISIE_ML4 <- function( verbose = 0, tolint = c(1E-16, 1E-10), island_ontogeny = NA, - jitter = 0) { + jitter = 0, + num_cycles = 1) { out2err <- data.frame( lambda_c = NA, @@ -194,7 +195,8 @@ DAISIE_ML4 <- function( CS_version = CS_version, abstolint = tolint[1], reltolint = tolint[2], - jitter = jitter + jitter = jitter, + num_cycles = num_cycles ) if (out$conv != 0) { cat( diff --git a/R/DAISIE_ML_CS.R b/R/DAISIE_ML_CS.R index 1384c529..350cdac2 100644 --- a/R/DAISIE_ML_CS.R +++ b/R/DAISIE_ML_CS.R @@ -182,7 +182,8 @@ DAISIE_ML_CS <- DAISIE_ML <- function( CS_version = 1, verbose = 0, tolint = c(1E-16, 1E-10), - jitter = 0) { + jitter = 0, + num_cycles = 1) { if (datatype == "single") { if (is.na(island_ontogeny)) { @@ -203,7 +204,8 @@ DAISIE_ML_CS <- DAISIE_ML <- function( CS_version = CS_version, verbose = verbose, tolint = tolint, - jitter = jitter) + jitter = jitter, + num_cycles = num_cycles) } else { out <- DAISIE_ML1(datalist = datalist, @@ -226,7 +228,8 @@ DAISIE_ML_CS <- DAISIE_ML <- function( CS_version = CS_version, verbose = verbose, tolint = tolint, - jitter = jitter) + jitter = jitter, + num_cycles = num_cycles) } } else { @@ -246,7 +249,8 @@ DAISIE_ML_CS <- DAISIE_ML <- function( CS_version = CS_version, verbose = verbose, tolint = tolint, - jitter = jitter) + jitter = jitter, + num_cycles = num_cycles) } } else { @@ -265,7 +269,8 @@ DAISIE_ML_CS <- DAISIE_ML <- function( optimmethod = optimmethod, verbose = verbose, tolint = tolint, - jitter = jitter) + jitter = jitter, + num_cycles = num_cycles) } return(out) } diff --git a/R/DAISIE_ML_IW.R b/R/DAISIE_ML_IW.R index 234cabd1..53738f96 100644 --- a/R/DAISIE_ML_IW.R +++ b/R/DAISIE_ML_IW.R @@ -90,12 +90,13 @@ DAISIE_ML_IW <- function( optimmethod = "subplex", verbose = 0, tolint = c(1E-16, 1E-14), - jitter = 0) { + jitter = 0, + num_cycles = 1) { options(warn = -1) out2err <- data.frame(lambda_c = NA, mu = NA, K = NA, gamma = NA, lambda_a = NA, loglik = NA, df = NA, conv = NA) out2err <- invisible(out2err) if (is.null(datalist[[1]]$brts_table)) { - datalist <- Add_brt_table(datalist) + datalist <- add_brt_table(datalist) } np <- datalist[[1]]$not_present if (is.null(np)) { @@ -153,7 +154,8 @@ DAISIE_ML_IW <- function( methode = methode, abstolint = tolint[1], reltolint = tolint[2], - jitter = jitter + jitter = jitter, + num_cycles = num_cycles ) if (out$conv != 0) { cat("Optimization has not converged. Try again with different initial values.\n") diff --git a/R/DAISIE_MW_ML.R b/R/DAISIE_MW_ML.R index 2a44832c..3b2ae608 100644 --- a/R/DAISIE_MW_ML.R +++ b/R/DAISIE_MW_ML.R @@ -218,6 +218,7 @@ DAISIE_MW_loglik_choosepar = function( #' the 11th parameter is the d0 for colonization and the 12th is the d0 for #' anagenesis. #' +#' @inheritParams default_params_doc #' @param datalist Data object containing information on colonisation and #' branching times of species for several islands or archipelagos, as well as the area, #' isolation and age of each of the islands/archipelagos. See data(archipelagos41) for @@ -267,8 +268,6 @@ DAISIE_MW_loglik_choosepar = function( #' @param optimmethod Method used in likelihood optimization. Default is #' "subplex" (see subplex package). Alternative is 'simplex' which was the #' method in previous versions. -#' @param CS_version For internal testing purposes only. Default is 1, the -#' original DAISIE code. #' @param verbose sets whether parameters and likelihood should be printed (1) #' or not (0) #' @param tolint Vector of two elements containing the absolute and relative @@ -304,6 +303,8 @@ DAISIE_MW_loglik_choosepar = function( #' machine. #' @param cpus Number of cpus used in parallel computing. Default is 3. Will #' not have an effect if parallel = 'no'. +#' @param num_cycles The number of cycles the optimizer will go through. +#' Default is 1. #' @return The output is a dataframe containing estimated parameters and #' maximum loglikelihood. #' \item{lambda_c0}{ gives the maximum likelihood estimate of lambda^c, @@ -385,7 +386,8 @@ DAISIE_MW_ML = function( distance_type = 'continent', distance_dep = 'power', parallel = 'local', - cpus = 3 + cpus = 3, + num_cycles = 1 ) { options(warn=-1) @@ -443,7 +445,24 @@ DAISIE_MW_ML = function( } cat("Optimizing the likelihood - this may take a while.","\n") utils::flush.console() - out = DDD::optimizer(optimmethod = optimmethod,optimpars = optimpars,fun = DAISIE_MW_loglik_choosepar,trparsopt = trparsopt,idparsopt = idparsopt,trparsfix = trparsfix,idparsfix = idparsfix,pars2 = pars2,datalist = datalist,methode = methode,CS_version = CS_version,abstolint = tolint[1],reltolint = tolint[2],distance_type = distance_type,parallel = parallel,cpus = cpus,distance_dep = distance_dep) + out = DDD::optimizer(optimmethod = optimmethod, + optimpars = optimpars, + fun = DAISIE_MW_loglik_choosepar, + trparsopt = trparsopt, + idparsopt = idparsopt, + trparsfix = trparsfix, + idparsfix = idparsfix, + pars2 = pars2, + datalist = datalist, + methode = methode, + CS_version = CS_version, + abstolint = tolint[1], + reltolint = tolint[2], + distance_type = distance_type, + parallel = parallel, + cpus = cpus, + distance_dep = distance_dep, + num_cycles = num_cycles) if(out$conv != 0) { cat("Optimization has not converged. Try again with different initial values.\n") diff --git a/R/DAISIE_SR_ML_CS.R b/R/DAISIE_SR_ML_CS.R index 74c5e5f2..d639e01e 100644 --- a/R/DAISIE_SR_ML_CS.R +++ b/R/DAISIE_SR_ML_CS.R @@ -64,6 +64,7 @@ DAISIE_SR_loglik_all_choosepar <- function( #' used, otherwise the information in the data is overruled. #' #' @aliases DAISIE_SR_ML_CS DAISIE_SR_ML +#' @inheritParams default_params_doc #' @param datalist Data object containing information on colonisation and #' branching times. This object can be generated using the DAISIE_dataprep #' function, which converts a user-specified data table into a data object, but @@ -101,34 +102,33 @@ DAISIE_SR_loglik_all_choosepar <- function( #' e.g. c(1,3) if lambda^c and K should not be optimized. #' @param parsfix The values of the parameters that should not be optimized #' @param idparsnoshift The ids of the parameters that should not be different -#' before and after the shift +#' before and after the shift. #' @param res Sets the maximum number of species for which a probability must -#' be computed, must be larger than the size of the largest clade +#' be computed, must be larger than the size of the largest clade #' @param ddmodel Sets the model of diversity-dependence: \cr \cr ddmodel = 0 : -#' no diversity dependence \cr ddmodel = 1 : linear dependence in speciation -#' rate \cr ddmodel = 11: linear dependence in speciation rate and in -#' immigration rate \cr ddmodel = 2 : exponential dependence in speciation -#' rate\cr ddmodel = 21: exponential dependence in speciation rate and in -#' immigration rate\cr +#' no diversity dependence \cr ddmodel = 1 : linear dependence in speciation +#' rate \cr ddmodel = 11: linear dependence in speciation rate and in +#' immigration rate \cr ddmodel = 2 : exponential dependence in speciation +#' rate\cr ddmodel = 21: exponential dependence in speciation rate and in +#' immigration rate.\cr #' @param cond cond = 0 : conditioning on island age \cr cond = 1 : -#' conditioning on island age and non-extinction of the island biota \cr -#' @param island_ontogeny type of island ontonogeny. If NA, then constant ontogeny is assumed +#' conditioning on island age and non-extinction of the island biota \cr +#' @param island_ontogeny type of island ontonogeny. If NA, then constant +#' ontogeny is assumed. #' @param tol Sets the tolerances in the optimization. Consists of: \cr reltolx -#' = relative tolerance of parameter values in optimization \cr reltolf = -#' relative tolerance of function value in optimization \cr abstolx = absolute -#' tolerance of parameter values in optimization -#' @param maxiter Sets the maximum number of iterations in the optimization +#' = relative tolerance of parameter values in optimization \cr reltolf = +#' relative tolerance of function value in optimization \cr abstolx = absolute +#' tolerance of parameter values in optimization. +#' @param maxiter Sets the maximum number of iterations in the optimization. #' @param methode Method of the ODE-solver. See package deSolve for details. -#' Default is "lsodes" +#' Default is "lsodes" #' @param optimmethod Method used in likelihood optimization. Default is -#' "subplex" (see subplex package). Alternative is 'simplex' which was the -#' method in previous versions. -#' @param CS_version For internal testing purposes only. Default is 1, the -#' original DAISIE code. +#' "subplex" (see subplex package). Alternative is 'simplex' which was the +#' method in previous versions. #' @param verbose sets whether parameters and likelihood should be printed (1) -#' or not (0) +#' or not (0). #' @param tolint Vector of two elements containing the absolute and relative -#' tolerance of the integration +#' tolerance of the integration. #' @param jitter Numeric for \code{\link[DDD]{optimizer}()}. Jitters the #' parameters being optimized by the specified amount which should be very #' small, e.g. 1e-5. Jitter when \code{link[subplex]{subplex}()} produces @@ -285,7 +285,8 @@ DAISIE_SR_ML_CS <- DAISIE_SR_ML <- function( CS_version = 1, verbose = 0, tolint = c(1E-16, 1E-10), - jitter = 0) { + jitter = 0, + num_cycles = 1) { # datalist = list of all data: branching times, status of clade, and numnber of missing species # datalist[[,]][1] = list of branching times (positive, from present to past) # - max(brts) = age of the island @@ -399,7 +400,8 @@ DAISIE_SR_ML_CS <- DAISIE_SR_ML <- function( CS_version = CS_version, abstolint = tolint[1], reltolint = tolint[2], - jitter = jitter + jitter = jitter, + num_cycles = num_cycles ) if (out$conv != 0) { cat("Optimization has not converged. Try again with different initial values.\n") diff --git a/R/DAISIE_SR_loglik_CS.R b/R/DAISIE_SR_loglik_CS.R index ffde67bf..9e142006 100644 --- a/R/DAISIE_SR_loglik_CS.R +++ b/R/DAISIE_SR_loglik_CS.R @@ -180,7 +180,8 @@ DAISIE_SR_loglik_CS_M1 <- DAISIE_SR_loglik <- function( } if (stac == 2 || stac == 3 || stac == 4) { gamvec <- divdepvecproc(pars1, lx, k1, ddep * (ddep == 11 | ddep == 21), brts[2], "col") - probs[(2 * lx + 1):(3 * lx)] <- gamvec[1:lx] * probs[1:lx] + probs[(2 * lx + 1):(3 * lx)] <- gamvec[1:lx] * probs[1:lx] + + gamvec[2:(lx + 1)] * probs[(lx + 1):(2 * lx)] probs[1:(2 * lx)] <- 0 k1 <- 1 y <- odeproc(probs, brts[2:3], DAISIE_loglik_rhs2, c(pars1, k1, ddep), rtol = reltolint, atol = abstolint, method = methode) @@ -257,6 +258,7 @@ DAISIE_SR_loglik_CS_M1 <- DAISIE_SR_loglik <- function( #' #' @aliases DAISIE_SR_loglik_CS DAISIE_SR_loglik_all #' +#' @inheritParams default_params_doc #' @param pars1 Contains the model parameters: \cr \cr #' \code{pars1[1]} #' corresponds to lambda^c (cladogenesis rate) \cr @@ -324,8 +326,6 @@ DAISIE_SR_loglik_CS_M1 <- DAISIE_SR_loglik <- function( #' applicable for endemic clades) \cr #' @param methode Method of the ODE-solver. See package deSolve for details. #' Default is "lsodes" -#' @param CS_version For internal testing purposes only. Default is 1, the -#' original DAISIE code. #' @param abstolint Absolute tolerance of the integration #' @param verbose Logical controling if progress is printed to console. #' @param reltolint Relative tolerance of the integration diff --git a/R/DAISIE_dataprep.R b/R/DAISIE_dataprep.R index 8c199638..fd5823d9 100644 --- a/R/DAISIE_dataprep.R +++ b/R/DAISIE_dataprep.R @@ -335,8 +335,7 @@ DAISIE_dataprep = function(datatable, if (length(which(is.na(unlist(datalist)[which(names(unlist(datalist)) == 'stac')]) == TRUE)) > 0) { - stop(paste("The status of one or more lineages is incorrectly spelled in - the source table and has not been assigned.")) + stop("The status of one or more lineages is incorrectly spelled in the source table and has not been assigned.") } return(datalist) diff --git a/R/DAISIE_format_IW.R b/R/DAISIE_format_IW.R index f72f9742..2ed32530 100644 --- a/R/DAISIE_format_IW.R +++ b/R/DAISIE_format_IW.R @@ -77,16 +77,20 @@ DAISIE_format_IW <- function(island_replicates, stt_all = stt_all ) } else { + taxon_list_size <- length(the_island$taxon_list) island_list[[1]] <- list( island_age = totaltime, - not_present = M - length(the_island$taxon_list), + not_present = M - taxon_list_size, stt_all = stt_all ) - for (y in 1:length(the_island$taxon_list)) { - island_list[[y + 1]] <- the_island$taxon_list[[y]] + if (taxon_list_size != 0) { + for (y in seq_len(taxon_list_size)) { + island_list[[y + 1]] <- the_island$taxon_list[[y]] + } } } - island_list <- Add_brt_table(island_list) + + island_list <- add_brt_table(island_list) several_islands[[rep]] <- island_list if (verbose == TRUE) { @@ -153,7 +157,7 @@ DAISIE_format_IW_trait <- function(island_replicates, } } - island_list = Add_brt_table(island_list) + island_list <- add_brt_table(island_list) several_islands[[rep]] = island_list if (verbose) { @@ -169,12 +173,12 @@ DAISIE_format_IW_trait <- function(island_replicates, return(several_islands) } -Add_brt_table <- function(island) { +add_brt_table <- function(island, full_table = FALSE) { island_age <- island[[1]]$island_age island_top <- island[[1]] if (length(island) == 1) { - brts_table <- matrix(ncol = 4, nrow = 1) - brts_table[1, ] <- c(island_age, 0, 0, NA) + brts_table <- matrix(ncol = 5, nrow = 1) + brts_table[1, ] <- c(island_age, 0, 0, NA, NA) island[[1]]$brts_table <- brts_table } else { island_top <- island[[1]] @@ -190,8 +194,8 @@ Add_brt_table <- function(island) { stac1_5s <- sort(c(stac1s, stac5s)) if (length(stac1_5s) != 0) { if (length(stac1_5s) == length(island)) { - brts_table <- matrix(ncol = 4, nrow = 1) - brts_table[1, ] <- c(island_age, 0, 0, NA) + brts_table <- matrix(ncol = 5, nrow = 1) + brts_table[1, ] <- c(island_age, 0, 0, NA, NA) island_no_stac1or5 <- NULL } else { island_no_stac1or5 <- island[-stac1_5s] @@ -206,21 +210,41 @@ Add_brt_table <- function(island) { btimes[[i]] <- island_no_stac1or5[[i]]$branching_times[-1] } brts <- rev(sort(unlist(btimes))) - brts_IWK <- matrix(ncol = 4, nrow = length(brts)) + brts_IWK <- NULL pos1 <- 0 + j <- 1 for (i in 1:length(btimes)) { the_brts <- btimes[[i]] the_stac <- island_no_stac1or5[[i]]$stac pos2 <- pos1 + length(the_brts) - brts_IWK[(pos1 + 1):pos2, 1] <- the_brts - brts_IWK[(pos1 + 1):pos2, 2] <- i - brts_IWK[(pos1 + 1):pos2, 3] <- seq(1, length(the_brts)) - brts_IWK[(pos1 + 1):pos2, 4] <- (the_stac == 2) + + ff <- matrix(ncol = 5, nrow = pos2 - pos1) + ff[1:(pos2 - pos1), 1] <- the_brts + ff[1:(pos2 - pos1), 2] <- i + ff[1:(pos2 - pos1), 3] <- seq(1, length(the_brts)) + ff[1:(pos2 - pos1), 4] <- (the_stac == 2) + (the_stac == 3) + (the_stac == 4) * 0 + ff[1:(pos2 - pos1), 5] <- NA + brts_IWK <- rbind(brts_IWK,ff) pos1 <- pos2 + j <- j + 1 + if( !is.null(island[[i]]$all_colonisations) & full_table == TRUE) { + for (k in 1:length(island[[i]]$all_colonisations)) { + the_brts <- island[[i]]$all_colonisations[[k]]$event_times[-1] + pos2 <- pos1 + length(the_brts) + ff <- matrix(ncol = 5, nrow = pos2 - pos1 + 1) + ff[1:(pos2 - pos1), 1] <- the_brts + ff[1:(pos2 - pos1), 2] <- j + ff[1:(pos2 - pos1), 3] <- seq(1, length(the_brts)) + ff[1:(pos2 - pos1), 4] <- NA + ff[1:(pos2 - pos1), 5] <- j - 1 + brts_IWK <- rbind(brts_IWK,ff) + pos1 <- pos2 + j <- j + 1 + } + } } brts_table <- brts_IWK[rev(order(brts_IWK[, 1])), ] - brts_table <- rbind(c(island_age, 0, 0, NA), brts_table) + brts_table <- rbind(c(island_age, 0, 0, NA, NA), brts_table) } island_top$brts_table <- brts_table if (length(stac1_5s) != 0) { @@ -232,6 +256,6 @@ Add_brt_table <- function(island) { } island <- append(list(island_top), island) } - colnames(island[[1]]$brts_table) <- c("brt", "clade", "event", "endemic") + colnames(island[[1]]$brts_table) <- c("brt", "clade", "event", "endemic", "col") return(island) } diff --git a/R/DAISIE_loglik_CS.R b/R/DAISIE_loglik_CS.R index 831a07f4..792d3543 100644 --- a/R/DAISIE_loglik_CS.R +++ b/R/DAISIE_loglik_CS.R @@ -82,15 +82,19 @@ DAISIE_loglik_rhs <- function(t, x, parsvec) { ix3 = nil2lx ix4 = nil2lx-2 - dx1 = laavec[il1 + 1] * xx2[ix1] + lacvec[il4 + 1] * xx2[ix4] + muvec[il2 + 1] * xx2[ix3] + - lacvec[il1] * nn[in1] * xx1[ix1] + muvec[il2] * nn[in2] * xx1[ix2] + + dx1 = laavec[il1 + 1] * xx2[ix1] + + lacvec[il4 + 1] * xx2[ix4] + + muvec[il2 + 1] * xx2[ix3] + + lacvec[il1] * nn[in1] * xx1[ix1] + + muvec[il2] * nn[in2] * xx1[ix2] + -(muvec[il3] + lacvec[il3]) * nn[in3] * xx1[ix3] + -gamvec[il3] * xx1[ix3] dx1[1] = dx1[1] + laavec[il3[1]] * xx3 * (kk == 1) dx1[2] = dx1[2] + 2 * lacvec[il3[1]] * xx3 * (kk == 1) dx2 = gamvec[il3] * xx1[ix3] + - lacvec[il1 + 1] * nn[in1] * xx2[ix1] + muvec[il2 + 1] * nn[in2] * xx2[ix2] + + lacvec[il1 + 1] * nn[in1] * xx2[ix1] + + muvec[il2 + 1] * nn[in2] * xx2[ix2] + -(muvec[il3 + 1] + lacvec[il3 + 1]) * nn[in3 + 1] * xx2[ix3] + -laavec[il3 + 1] * xx2[ix3] @@ -144,9 +148,13 @@ DAISIE_loglik_rhs2 <- function(t, x, parsvec) { # outflow: # all events with n+k species present dx1 = (laavec[il3] * xx3[ix3] + 2 * lacvec[il1] * xx3[ix1]) * (kk == 1) + - laavec[il1 + 1] * xx2[ix1] + lacvec[il4 + 1] * xx2[ix4] + muvec[il2 + 1] * xx2[ix3] + - lacvec[il1] * nn[in1] * xx1[ix1] + muvec[il2] * nn[in2] * xx1[ix2] + - -(muvec[il3] + lacvec[il3]) * nn[in3] * xx1[ix3] - gamvec[il3] * xx1[ix3] + laavec[il1 + 1] * xx2[ix1] + + lacvec[il4 + 1] * xx2[ix4] + + muvec[il2 + 1] * xx2[ix3] + + lacvec[il1] * nn[in1] * xx1[ix1] + + muvec[il2] * nn[in2] * xx1[ix2] + + -(muvec[il3] + lacvec[il3]) * nn[in3] * xx1[ix3] + + -gamvec[il3] * xx1[ix3] # inflow: # immigration when there are n+k species: Q^k,n -> Q^M,k_n; @@ -157,7 +165,8 @@ DAISIE_loglik_rhs2 <- function(t, x, parsvec) { # outflow: # all events with n+k+1 species present dx2 <- gamvec[il3] * xx1[ix3] + - lacvec[il1 + 1] * nn[in1] * xx2[ix1] + muvec[il2 + 1] * nn[in2] * xx2[ix2] + + lacvec[il1 + 1] * nn[in1] * xx2[ix1] + + muvec[il2 + 1] * nn[in2] * xx2[ix2] + -(muvec[il3 + 1] + lacvec[il3 + 1]) * nn[in3 + 1] * xx2[ix3] + -laavec[il3 + 1] * xx2[ix3] @@ -169,7 +178,8 @@ DAISIE_loglik_rhs2 <- function(t, x, parsvec) { # n+k+1 species present # outflow: # all events with n+k species present - dx3 <- lacvec[il1] * nn[in4] * xx3[ix1] + muvec[il2] * nn[in2] * xx3[ix2] + + dx3 <- lacvec[il1] * nn[in4] * xx3[ix1] + + muvec[il2] * nn[in2] * xx3[ix2] + -(lacvec[il3] + muvec[il3]) * nn[in3] * xx3[ix3] + -(laavec[il3] + gamvec[il3]) * xx3[ix3] @@ -441,8 +451,9 @@ DAISIE_loglik_CS_M1 <- DAISIE_loglik <- function(pars1, } if (stac == 2 || stac == 3 || stac == 4) { t <- brts[2] - gamvec = divdepvec(gam,c(pars1_in_divdepvec_call,t,0),lx,k1,ddep * (ddep == 11 | ddep == 21),island_ontogeny) # Problem may be here 30/3 - probs[(2 * lx + 1):(3 * lx)] = gamvec[1:lx] * probs[1:lx] + gamvec = divdepvec(gam,c(pars1_in_divdepvec_call,t,0),lx,k1,ddep * (ddep == 11 | ddep == 21),island_ontogeny) + probs[(2 * lx + 1):(3 * lx)] = gamvec[1:lx] * probs[1:lx] + + gamvec[2:(lx + 1)] * probs[(lx + 1):(2 * lx)] probs[1:(2 * lx)] = 0 k1 = 1 #y = deSolve::ode(probs,c(brts[2:3]),DAISIE_loglik_rhs2,c(pars1,k1,ddep),rtol = reltolint,atol = abstolint,method = methode) @@ -548,9 +559,10 @@ DAISIE_loglik_CS_M1 <- DAISIE_loglik <- function(pars1, return(loglik) } -DAISIE_loglik_CS_choice = function( +DAISIE_loglik_CS_choice <- function( pars1, pars2, + datalist = NULL, brts, stac, missnumspec, @@ -590,6 +602,7 @@ DAISIE_loglik_CS_choice = function( loglik <- DAISIE_loglik_IW_M1( pars1 = pars1, pars2 = pars2, + datalist = datalist, brts = brts, stac = stac, missnumspec = missnumspec, @@ -610,6 +623,7 @@ DAISIE_loglik_CS_choice = function( #' @description Computes the loglikelihood of the DAISIE model with clade-specific #' diversity-dependence given colonization and branching times for lineages on #' an island, and a set of model parameters. The output is a loglikelihood value +#' @inheritParams default_params_doc #' @param pars1 Contains the model parameters: \cr \cr #' \code{pars1[1]} corresponds to lambda^c (cladogenesis rate) \cr #' \code{pars1[2]} corresponds to mu (extinction rate) \cr @@ -675,13 +689,13 @@ DAISIE_loglik_CS_choice = function( #' * Endemic_Singleton_MaxAge: 5 \cr #' * Endemic_Clade_MaxAge: 6 \cr #' * Endemic&Non_Endemic_Clade_MaxAge: 7 \cr \cr +#' * Non_endemic_MaxAge_MinAge: 8 \cr +#' * Endemic_Singleton_MaxAge_MinAge: 9 \cr #' \code{$missing_species} - number of island species that were not sampled for #' particular clade (only applicable for endemic clades) \cr #' \code{$type1or2} - whether the colonist belongs to type 1 or type 2 \cr #' @param methode Method of the ODE-solver. See package deSolve for details. #' Default is "lsodes" -#' @param CS_version For internal testing purposes only. Default is 1, the -#' original DAISIE code. #' @param abstolint Absolute tolerance of the integration #' @param reltolint Relative tolerance of the integration #' @return The loglikelihood @@ -710,8 +724,8 @@ DAISIE_loglik_CS <- DAISIE_loglik_all <- function(pars1, CS_version = 1, abstolint = 1E-16, reltolint = 1E-10) { - pars1 = as.numeric(pars1) - cond = pars2[3] + pars1 <- as.numeric(pars1) + cond <- pars2[3] endpars1 <- 5 if(length(pars1) == 5 | !is.na(pars2[5])) # Normal no ont case @@ -806,23 +820,30 @@ DAISIE_loglik_CS <- DAISIE_loglik_all <- function(pars1, { if(datalist[[i]]$type1or2 == 1) { - pars = pars1[1:endpars1] + pars <- pars1[1:endpars1] } else { pars <- pars1[6:10] } loglik <- loglik + DAISIE_loglik_CS_choice( pars1 = pars, pars2 = pars2, + datalist = datalist[[i]], brts = datalist[[i]]$branching_times, stac = datalist[[i]]$stac, missnumspec = datalist[[i]]$missing_species, methode = methode, CS_version = CS_version, abstolint = abstolint, - reltolint = reltolint - ) + reltolint = reltolint) } } + if (pars2[4] >= 1) { + s1 <- sprintf("Parameters: ") + s2 <- sprintf("%f ",pars) + s3 <- sprintf(", Loglikelihood: %f", loglik) + cat(s1, s2, s3, "\n", sep = "") + utils::flush.console() + } return(loglik) } @@ -922,9 +943,9 @@ logcondprob <- function(numcolmin, numimm, logp0, fac = 2) { logcond <- 0 if(numcolmin >= 1) { if(numcolmin == 1 && length(logp0) == 2) { - cat('With two types, conditioning on at least one colonization - implies at least two colonizations. Therefore, the minimum - number of colonizations is changed to 2.\n') + message('With two types, conditioning on at least one colonization + implies at least two colonizations. Therefore, the minimum + number of colonizations is changed to 2.\n') numcolmin <- 2 } lognotp0 <- log1p(-exp(logp0)) @@ -936,23 +957,21 @@ logcondprob <- function(numcolmin, numimm, logp0, fac = 2) { pc <- exp(logpc) if(length(logp0) == 2) { condprob <- pc[1,1] + pc[1,2] - pc[1,1] * pc[1,2] - #condprob <- sum(pc[1,1] * pc[,2]) + sum(pc[1,2] * pc[,1]) - pc[1,1] * pc[1,2] if(numcolmin > 2) { for(i in 2:(numcolmin - 1)) { - condprob <- condprob + sum(pc[2:i,1] * pc[i:2,2]) + condprob <- condprob + sum(pc[2:i,1] * pc[i:2,2]) } } if(condprob >= 1) { logcond <- log(sum(pc[2:numcolmin,1] * pc[numcolmin:2,2])) - cat('A simple approximation of logcond must be made. Results may be unreliable.\n') + message('A simple approximation of logcond must be made. Results may be unreliable.\n') } else { logcond <- log1p(-condprob) } } else { - #if(sum(pc[-(numcolmin + 1)]) >= 1) { if(sum(pc) >= 1) { logcond <- log(sum(pc[(numcolmin + 1):(maxi + 1)])) - cat('An approximation of logcond must be made. Results may be unreliable.\n') + message('An approximation of logcond must be made. Results may be unreliable.\n') } else { logcond <- log1p(-sum(pc[-((numcolmin + 1):(maxi + 1))])) } diff --git a/R/DAISIE_loglik_IW.R b/R/DAISIE_loglik_IW.R index 4f82ce8c..19ca5993 100644 --- a/R/DAISIE_loglik_IW.R +++ b/R/DAISIE_loglik_IW.R @@ -41,48 +41,46 @@ kimat <- function(dec2binmatk) { return(ki) } -kmini <- function(dec2binmatk, lxm1, lxm2, lxe, sysdim = dim(dec2binmatk)[1]) { +create_l0ki <- function(dec2binmatk, lxm, lxe, sysdim = dim(dec2binmatk)[1]) { posc <- Matrix::rowSums(dec2binmatk) negc <- log2(sysdim) - posc - kmin <- rep(negc, each = lxm1 * lxm2 * lxe) - dim(kmin) <- c(lxm1, lxm2, lxe, sysdim) + l0 <- rep(negc, each = lxm * lxe) + dim(l0) <- c(lxm, lxe, sysdim) ki <- kimat(dec2binmatk) - res <- list(kmin = kmin, ki = ki) + res <- list(l0 = l0, ki = ki) return(res) } -nndivdep <- function(lxm1, lxm2, lxe, sysdim, Kprime, M, k, l) { - nnl <- c(0, 0:(lxm1 + 1)) - nnm <- c(0, 0:(lxm2 + 1)) +nndivdep <- function(lxm, lxe, sysdim, Kprime, k, M, l0) { + nnm <- c(0, 0:(lxm + 1)) nne <- c(0, 0, 0:(lxe + 1)) - lnnl <- length(nnl) lnnm <- length(nnm) lnne <- length(nne) - nil2lxm1 <- 2:(lxm1 + 1) - nil2lxm2 <- 2:(lxm2 + 1) + nil2lxm <- 2:(lxm + 1) nil2lxe <- 3:(lxe + 2) - nn <- rowSums(expand.grid(n1 = nnl, n2 = nnm, n3 = nne)) - dim(nn) <- c(lnnl, lnnm, lnne) + nn <- rowSums(expand.grid(n1 = nnm, n2 = nne)) + dim(nn) <- c(lnnm, lnne) nn <- replicate(sysdim, nn) - nilm1 <- rep(1, lxm1) - nilm2 <- rep(1, lxm2) + nilm <- rep(1, lxm) nile <- rep(1, lxe) allc <- 1:sysdim - divdepfac <- pmax(array(0, dim = c(lxm1 + 3, lxm2 + 3, lxe + 4, sysdim)), + divdepfac <- pmax(array(0, dim = c(lxm + 3, lxe + 4, sysdim)), 1 - (nn + k) / Kprime) - divdepfacmin1 <- pmax(array(0, dim = c(lxm1 + 3, lxm2 + 3, lxe + 4, sysdim)), + divdepfacmin1 <- pmax(array(0, dim = c(lxm + 3, lxe + 4, sysdim)), 1 - (nn + k - 1) / Kprime) - divdepfac <- divdepfac[nil2lxm1, nil2lxm2, nil2lxe, allc] - divdepfacmin1 <- divdepfacmin1[nil2lxm1, nil2lxm2, nil2lxe, allc] - Mminm <- M - nn[nil2lxm1, nil2lxm2, nile, allc] - lminm1 <- l - nn[nil2lxm1, nilm2, nile, allc] - Mminlminm2 <- M - l - nn[nilm1, nil2lxm2, nile, allc] + divdepfacplus1 <- pmax(array(0, dim = c(lxm + 3, lxe + 4, sysdim)), + 1 - (nn + k + 1) / Kprime) + divdepfac <- divdepfac[nil2lxm, nil2lxe, allc] + divdepfacmin1 <- divdepfacmin1[nil2lxm, nil2lxe, allc] + divdepfacplus1 <- divdepfacplus1[nil2lxm, nil2lxe, allc] + mfac <- (nn[nil2lxm,nile,allc] + 1)/(M - l0) + oneminmfac <- (M - nn[nil2lxm,nile,allc] - l0)/(M - l0) res <- list(nn = nn, divdepfac = divdepfac, divdepfacmin1 = divdepfacmin1, - Mminm = Mminm, - lminm1 = lminm1, - Mminlminm2 = Mminlminm2) + divdepfacplus1 = divdepfacplus1, + mfac = mfac, + oneminmfac = oneminmfac) return(res) } @@ -96,151 +94,74 @@ selectrows <- function(sysdim, order) { return(mat) } - -DAISIE_odeint_iw_pars <- function(pars) -{ - lac = pars[[1]][1] - mu = pars[[1]][2] - Kprime = pars[[1]][3] - gam = pars[[1]][4] - laa = pars[[1]][5] - M = pars[[1]][6] - k = pars[[2]] - ddep = pars[[3]] - lxm1 = pars[[4]]$lxm1 - lxm2 = pars[[4]]$lxm2 - lxe = pars[[4]]$lxe - sysdim = pars[[4]]$sysdim - kmin = pars[[5]]$kmin - kplus = k - kmin - ki = pars[[5]]$ki - nn = pars[[6]]$nn - divdepfac = pars[[6]]$divdepfac - divdepfacmin1 = pars[[6]]$divdepfacmin1 - Mminm = pars[[6]]$Mminm - Mminm[Mminm < 0] = 0 - lminm1minkminplus1 = pars[[6]]$lminm1 - kmin + 1 - lminm1minkminplus1[lminm1minkminplus1 < 0] = 0 - Mminlminm2plus1 = pars[[6]]$Mminlminm2 + 1 - Mminlminm2plus1[Mminlminm2plus1 < 0] = 0 - - nil2lxm1 = 2:(lxm1 + 1) - nil2lxm2 = 2:(lxm2 + 1) - nil2lxe = 3:(lxe + 2) - nilm1 = rep(1,lxm1) - nile = rep(1,lxe) - nilm2 = rep(1,lxm2) - - allc = 1:sysdim - if(sysdim == 1 & lxm1 > 1) - { - dim(Mminm) = c(lxm1,lxm2,lxe) - } else - if(sysdim == 1 & lxm1 == 1) - { - dim(Mminm) = c(lxm2,lxe) - } else - if (sysdim == 1 & lxm1 == 1) { - dim(Mminm) <- c(lxm2, lxe) - } else - if (sysdim > 1 & lxm1 == 1) { - dim(Mminm) <- c(lxm2, lxe, allc) - } - - CP = list( - c1 = (gam * divdepfacmin1 * lminm1minkminplus1), - c2 = (gam * divdepfacmin1 * Mminlminm2plus1), - c3 = (mu * nn[nil2lxm1 + 1, nilm2, nile, allc]), - c4 = (mu * nn[nilm1, nil2lxm2 + 1, nile, allc]), - c5 = (mu * nn[nilm1, nilm2, nil2lxe + 1, allc]), - c6 = (lac * divdepfacmin1 * nn[nil2lxm1 + 1, nilm2, nile, allc]), - c7 = (lac * divdepfacmin1 * nn[nilm1, nil2lxm2 + 1, nile, allc]), - c89 = ((lac * divdepfacmin1 * nn[nilm1, nilm2, nil2lxe - 1, allc]) + (2 * kplus * lac * divdepfacmin1)), - c10 = (laa * nn[nil2lxm1 + 1, nilm2, nile, allc]), - c11 = (laa * nn[nilm1, nil2lxm2 + 1, nile, allc]), - c12 = (-(laa * (nn[nil2lxm1, nil2lxm2, nile, allc] + kmin) + (gam * divdepfac * Mminm) + - ((lac * divdepfac + mu) * (nn[nil2lxm1, nil2lxm2, nil2lxe, allc] + k)))), - c13 = (2 * lac * divdepfacmin1), - ki = ki, +DAISIE_IW_pars <- function(parslist) { + lac <- parslist$pars[1] + mu <- parslist$pars[2] + Kprime <- parslist$pars[3] + gam <- parslist$pars[4] + laa <- parslist$pars[5] + M <- parslist$pars[6] + k <- parslist$k + ddep <- parslist$ddep + lxm <- parslist$dime$lxm + lxe <- parslist$dime$lxe + sysdim <- parslist$dime$sysdim + l0 <- parslist$l0ki$l0 + ki <- parslist$l0ki$ki + nn <- parslist$nndd$nn + divdepfac <- parslist$nndd$divdepfac + divdepfacmin1 <- parslist$nndd$divdepfacmin1 + nil2lxm <- 2:(lxm + 1) + nil2lxe <- 3:(lxe + 2) + allc <- 1:sysdim + nilm <- rep(1, lxm) + nile <- rep(1, lxe) + cp <- list( laa = laa, - sysdim = sysdim + ki = ki, + lxm = lxm, + lxe = lxe, + sysdim = sysdim, + c1 = gam * divdepfacmin1 * pmax(0,M - l0 - nn[nil2lxm - 1, nile, allc]), + c2 = mu * nn[nil2lxm + 1, nile, allc], + c3 = mu * nn[nilm, nil2lxe + 1, allc], + c4 = laa * nn[nil2lxm + 1, nile, allc], + c5 = lac * divdepfacmin1 * nn[nil2lxm + 1, nile, allc], + c6 = lac * divdepfacmin1 * (2 * (k - l0) + nn[nilm, nil2lxe - 1, allc]), + c7 = gam * divdepfac * pmax(0,M - nn[nil2lxm, nile, allc]) + + mu * (k + nn[nil2lxm, nil2lxe, allc]) + + laa * (l0 + nn[nil2lxm, nile, allc]) + + lac * divdepfac * (k + nn[nil2lxm, nil2lxe, allc]), + c8 = 2 * lac * divdepfacmin1 ) - return(CP) + return(cp) } - -DAISIE_loglik_rhs_IW = function(t,x,pars) +DAISIE_loglik_rhs_IW <- function(t,x,cp) { - lac = pars[[1]][1] - mu = pars[[1]][2] - Kprime = pars[[1]][3] - gam = pars[[1]][4] - laa = pars[[1]][5] - M = pars[[1]][6] - k = pars[[2]] - ddep = pars[[3]] - lxm1 = pars[[4]]$lxm1 - lxm2 = pars[[4]]$lxm2 - lxe = pars[[4]]$lxe - sysdim = pars[[4]]$sysdim - kmin = pars[[5]]$kmin - kplus = k - kmin - ki = pars[[5]]$ki - nn = pars[[6]]$nn - divdepfac = pars[[6]]$divdepfac - divdepfacmin1 = pars[[6]]$divdepfacmin1 - Mminm = pars[[6]]$Mminm - Mminm[Mminm < 0] = 0 - lminm1minkminplus1 = pars[[6]]$lminm1 - kmin + 1 - lminm1minkminplus1[lminm1minkminplus1 < 0] = 0 - Mminlminm2plus1 = pars[[6]]$Mminlminm2 + 1 - Mminlminm2plus1[Mminlminm2plus1 < 0] = 0 - - dim(x) = c(lxm1,lxm2,lxe,sysdim) - xx = array(0,dim = c(lxm1+2,lxm2+2,lxe+3,sysdim)) - xx[2:(lxm1+1),2:(lxm2+1),3:(lxe+2),1:sysdim] = x - nil2lxm1 = 2:(lxm1 + 1) - nil2lxm2 = 2:(lxm2 + 1) - nil2lxe = 3:(lxe + 2) - allc = 1:sysdim - nilm1 = rep(1,lxm1) - nile = rep(1,lxe) - nilm2 = rep(1,lxm2) - - if(sysdim == 1 & lxm1 > 1) - { - dim(Mminm) = c(lxm1,lxm2,lxe) - } else - if(sysdim == 1 & lxm1 == 1) - { - dim(Mminm) = c(lxm2,lxe) - } else - if (sysdim == 1 & lxm1 == 1) { - dim(Mminm) <- c(lxm2, lxe) - } else - if (sysdim > 1 & lxm1 == 1) { - dim(Mminm) <- c(lxm2, lxe, allc) - } - - dx <- gam * divdepfacmin1 * lminm1minkminplus1 * xx[nil2lxm1 - 1, nil2lxm2, nil2lxe, allc] + #immigration - gam * divdepfacmin1 * Mminlminm2plus1 * xx[nil2lxm1, nil2lxm2 - 1, nil2lxe, allc] + #immigration - mu * nn[nil2lxm1 + 1, nilm2, nile, allc] * xx[nil2lxm1 + 1, nil2lxm2, nil2lxe, allc] + #extinction non-endemics - mu * nn[nilm1, nil2lxm2 + 1, nile, allc] * xx[nil2lxm1, nil2lxm2 + 1, nil2lxe, allc] + #extinction non-endemics - mu * nn[nilm1, nilm2, nil2lxe + 1, allc] * xx[nil2lxm1, nil2lxm2, nil2lxe + 1, allc] + #extinction endemics - lac * divdepfacmin1 * nn[nil2lxm1 + 1, nilm2, nile, allc] * xx[nil2lxm1 + 1, nil2lxm2, nil2lxe - 2, allc] + #cladogenesis non-endemics - lac * divdepfacmin1 * nn[nilm1, nil2lxm2 + 1, nile, allc] * xx[nil2lxm1, nil2lxm2 + 1, nil2lxe - 2, allc] + #cladogenesis non-endemics - lac * divdepfacmin1 * nn[nilm1, nilm2, nil2lxe - 1, allc] * xx[nil2lxm1, nil2lxm2, nil2lxe - 1, allc] + #cladogenesis endemics - 2 * kplus * lac * divdepfacmin1 * xx[nil2lxm1, nil2lxm2, nil2lxe - 1, allc] + #cladogenesis species in tree - laa * nn[nil2lxm1 + 1, nilm2, nile, allc] * xx[nil2lxm1 + 1, nil2lxm2, nil2lxe - 1, allc] + #anagenesis non-endemics - laa * nn[nilm1, nil2lxm2 + 1, nile, allc] * xx[nil2lxm1, nil2lxm2 + 1, nil2lxe - 1, allc] + #anagenesis non-endemics - -(laa * (nn[nil2lxm1, nil2lxm2, nile, allc] + kmin) + (gam * divdepfac * Mminm) + - (lac * divdepfac + mu) * (nn[nil2lxm1, nil2lxm2, nil2lxe, allc] + k)) * xx[nil2lxm1, nil2lxm2, nil2lxe, allc] + lxm <- cp$lxm + lxe <- cp$lxe + sysdim <- cp$sysdim + dim(x) <- c(lxm ,lxe, sysdim) + xx <- array(0,dim = c(lxm + 2, lxe + 3, sysdim)) + nil2lxm <- 2:(lxm + 1) + nil2lxe <- 3:(lxe + 2) + allc <- 1:sysdim + xx[nil2lxm, nil2lxe, allc] <- x + dx <- + cp$c1 * xx[nil2lxm - 1, nil2lxe, allc] + + cp$c2 * xx[nil2lxm + 1, nil2lxe, allc] + + cp$c3 * xx[nil2lxm, nil2lxe + 1, allc] + + cp$c4 * xx[nil2lxm + 1, nil2lxe - 1, allc] + + cp$c5 * xx[nil2lxm + 1, nil2lxe - 2, allc] + + cp$c6 * xx[nil2lxm, nil2lxe - 1, allc] - + cp$c7 * xx[nil2lxm, nil2lxe, allc] if (sysdim > 1) { dx <- dx + - laa * tensor::tensor(xx[nil2lxm1, nil2lxm2, nil2lxe, allc], ki, 4, 2) + # anagenesis in colonizing lineage - 2 * lac * divdepfacmin1 * tensor::tensor(xx[nil2lxm1, nil2lxm2, nil2lxe - 1, allc], ki, 4, 2) # cladogenesis in colonizing lineage + cp$laa * tensor::tensor(xx[nil2lxm, nil2lxe, allc], cp$ki, 3, 2) + + cp$c8 * tensor::tensor(xx[nil2lxm, nil2lxe - 1, allc], cp$ki, 3, 2) } - dim(dx) <- c(sysdim * lxm1 * lxm2 * lxe, 1) + dim(dx) <- c(lxm * lxe * sysdim, 1) return(list(dx)) } @@ -339,11 +260,12 @@ DAISIE_loglik_IW <- function( pars2[4] <- 0 } if (is.null(datalist[[1]]$brts_table)) { - datalist <- Add_brt_table(datalist) + datalist <- add_brt_table(datalist) } - brts <- c(-abs(datalist[[1]]$brts_table[,1]),0) - clade <- datalist[[1]]$brts_table[,2] - event <- datalist[[1]]$brts_table[,3] + brts <- c(-abs(datalist[[1]]$brts_table[,'brt']),0) + clade <- datalist[[1]]$brts_table[,'clade'] + event <- datalist[[1]]$brts_table[,'event'] + col <- datalist[[1]]$brts_table[,'col'] pars1 <- as.numeric(pars1) if(length(pars1) == 5) { @@ -359,8 +281,9 @@ DAISIE_loglik_IW <- function( return(loglik) } M <- length(datalist) - 1 + np + pars1[6] <- M } else - if (length(pars1) == 6) { + if(length(pars1) == 6) { M <- pars1[6] } else { cat("pars1 should contain 5 or 6 values.\n") @@ -383,7 +306,7 @@ DAISIE_loglik_IW <- function( } gam <- pars1[4] laa <- pars1[5] - l <- 0 + #l <- 0 if (min(pars1) < 0) { @@ -392,107 +315,136 @@ DAISIE_loglik_IW <- function( return(loglik) } - endemic = 0 - nonendemic1 = 0 - nonendemic2 = 0 - for(i in 2:length(datalist)) + if(length(datalist) > 1) for(i in 2:length(datalist)) { - endemic = endemic + (datalist[[i]]$stac == 5) - nonendemic1 = nonendemic1 + (datalist[[i]]$stac == 3) - nonendemic2 = nonendemic2 + (datalist[[i]]$stac == 1) + if(!datalist[[i]]$stac %in% c(0,2,4)) { + stop('IW does not work on data with unknown colonization times.') + } } - if((ddep == 1 | ddep == 11) & (ceiling(Kprime) < nonendemic1 + nonendemic2 + endemic + length(brts) - 2)) + if((ddep == 1 | ddep == 11) & (ceiling(Kprime) < length(brts) - 2)) { - if (verbose) { - cat('The proposed value of K is incompatible with the number of species - in the clade. Likelihood for this parameter set - will be set to -Inf. \n') - } + message('The proposed value of K is incompatible with the number of species + in the clade. Likelihood for this parameter set will be set to -Inf. \n') loglik <- -Inf return(loglik) } + if (ddep == 1 | ddep == 11) { lx <- min(1 + ceiling(Kprime), DDD::roundn(pars2[1]) ) } else { lx <- DDD::roundn(pars2[1]) } - lxm1 = max(clade) + 1 - lxm2 = min(lx,M + 1) - lxe = lx + lxm <- min(lx,M + 1) + lxe <- lx - sysdimchange = 1 - sysdim = 1 - totdim = lxm1 * lxm2 * lxe * sysdim - probs = rep(0,totdim) - probs[1] = 1 - loglik = 0 - expandvec = NULL - for(k in 0:(length(brts) - 2)) + if(M * (1 - exp((min(brts) * gam))) > 0.2 * lxm) { + message('With this colonization rate and system size setting, results may not be accurate.\n') + } + + sysdimchange <- 1 + sysdim <- 1 + totdim <- lxm * lxe * sysdim + probs <- rep(0, totdim) + probs[1] <- 1 + loglik <- 0 + expandvec <- NULL + for (k in 0:(length(brts) - 2)) { - if(pars2[4] == 2) + if (pars2[4] == 2) { cat(paste('k = ',k,', sysdim = ',sysdim,'\n',sep = '')) utils::flush.console() } - dime <- list(lxm1 = lxm1, lxm2 = lxm2, lxe = lxe, sysdim = sysdim) + dime <- list(lxm = lxm, lxe = lxe, sysdim = sysdim) if (sysdimchange == 1) { if (sysdim > 1) { dec2binmatk <- dec2binmat(log2(sysdim)) - kmi <- kmini(dec2binmatk, lxm1, lxm2, lxe, sysdim) + l0ki <- create_l0ki(dec2binmatk, lxm, lxe, sysdim) } else if (sysdim == 1) { - kmi <- list(kmin = 0, ki = NULL) + l0ki <- list(l0 = 0, ki = NULL) } sysdimchange <- 0 } - nndd = nndivdep(lxm1,lxm2,lxe,sysdim,Kprime,M,k,l) - parslist = list(pars = pars1,k = k,ddep = ddep,dime = dime,kmi = kmi,nndd = nndd) + nndd <- nndivdep(lxm = lxm, lxe = lxe, sysdim = sysdim, Kprime = Kprime, k = k, M = M, l0 = l0ki$l0) + parslist <- list(pars = pars1,k = k,ddep = ddep,dime = dime,l0ki = l0ki,nndd = nndd) + iw_parms = DAISIE_IW_pars(parslist) if (startsWith(methode, "odeint::")) { - probs = .Call("daisie_odeint_iw", probs, brts[(k + 1):(k + 2)], DAISIE_odeint_iw_pars(parslist), methode, abstolint, reltolint) - } - else { - y = deSolve::ode(y = probs,times = brts[(k + 1):(k + 2)],func = DAISIE_loglik_rhs_IW,parms = parslist,rtol = reltolint,atol = abstolint,method = methode) - probs = y[2,2:(totdim + 1)] + probs <- .Call("daisie_odeint_iw", probs, brts[(k + 1):(k + 2)], iw_parms, methode, abstolint, reltolint) + } else { + y <- deSolve::ode(y = probs, + times = brts[(k + 1):(k + 2)], + func = DAISIE_loglik_rhs_IW, + parms = iw_parms, + rtol = reltolint, + atol = abstolint, + method = methode) + probs <- y[2,2:(totdim + 1)] } - cp = checkprobs2(NA, loglik, probs, verbose); loglik = cp[[1]]; probs = cp[[2]] - dim(probs) = c(lxm1,lxm2,lxe,sysdim) + cp <- checkprobs2(NA, loglik, probs, verbose); loglik = cp[[1]]; probs = cp[[2]] + dim(probs) <- c(lxm, lxe, sysdim) if(k < (length(brts) - 2)) { - divdepfac = nndd$divdepfac - if(event[k + 2] == 1) + divdepfac <- nndd$divdepfac + divdepfacplus1 <- nndd$divdepfacplus1 + mfac <- nndd$mfac + oneminmfac <- nndd$oneminmfac + if(event[k + 2] == 1) #colonization { - Mminlminm2 = nndd$Mminlminm2 - Mminlminm2[Mminlminm2 < 0] = 0 - probs = gam * divdepfac * Mminlminm2 * probs[,,,1:sysdim] - probs = c(probs,rep(0,totdim)) - l = l + 1 - sysdim = sysdim * 2 - expandvec = c(expandvec,clade[k + 2]) - sysdimchange = 1 - } else + #l <- l + 1 + if(is.na(col[k + 2])) + { + test_for_colonization <- TRUE + } else + { + test_for_colonization <- (max(event[which(clade == col[k + 2])]) > 1) + } + if(test_for_colonization) # new colonization or recolonization after speciation + { + probs2 <- array(0,dim = dim(probs)) + probs2[1:(lxm - 1),,] <- probs[2:lxm,,] + #probs <- gam * divdepfac * probs[,,1:sysdim] + probs <- gam * divdepfac * oneminmfac * probs[,,1:sysdim] + + gam * divdepfacplus1 * mfac * probs2[,,1:sysdim] + probs <- c(probs,rep(0,totdim)) + sysdim <- sysdim * 2 + } else # recolonization without speciation + { + tocollapse <- which(expandvec == col[k + 2]) + sr <- selectrows(sysdim,2^(tocollapse - 1)) + probs[,,sr[,1]] <- 0 + probs <- gam * divdepfac * probs[,,1:sysdim] # + probs <- probs[,,sr[,2]] + expandvec <- expandvec[-tocollapse] + probs <- c(probs,rep(0,totdim/2)) + } + expandvec <- c(expandvec,clade[k + 2]) + sysdimchange <- 1 + } else # speciation { - probs = lac * divdepfac * probs[,,,1:sysdim] - if(event[k + 2] == 2) + probs <- lac * divdepfac * probs[,,1:sysdim] + if(event[k + 2] == 2) # first speciation in clade { - tocollapse = which(expandvec == clade[k + 2]) - sr = selectrows(sysdim,2^(tocollapse - 1)) - probs = probs[,,,sr[,1]] + probs[,,,sr[,2]] - sysdim = sysdim / 2 - dim(probs) = c(lxm1,lxm2,lxe,sysdim) - expandvec = expandvec[-tocollapse] - sysdimchange = 1 + tocollapse <- which(expandvec == clade[k + 2]) + sr <- selectrows(sysdim,2^(tocollapse - 1)) + probs <- probs[,,sr[,1]] + probs[,,sr[,2]] + sysdim <- sysdim / 2 + dim(probs) <- c(lxm,lxe,sysdim) + expandvec <- expandvec[-tocollapse] + sysdimchange <- 1 } } - cp = checkprobs2(NA, loglik, probs, verbose); loglik = cp[[1]]; probs = cp[[2]] - totdim = lxm1 * lxm2 * lxe * sysdim - dim(probs) = c(totdim,1) + cp <- checkprobs2(NA, loglik, probs, verbose); loglik <- cp[[1]]; probs <- cp[[2]] + totdim <- lxm * lxe * sysdim + dim(probs) <- c(totdim,1) } } - dim(probs) <- c(lxm1, lxm2, lxe, sysdim) + dim(probs) <- c(lxm, lxe, sysdim) expandedclades <- which(pracma::histc(clade, 1:length(clade))$cnt == 1) - status <- rep(0, lexpandedclades <- length(expandedclades)) + lexpandedclades <- length(expandedclades) + status <- rep(0, lexpandedclades) if (lexpandedclades > 0) { for (i in lexpandedclades:1) { if (datalist[[1 + expandedclades[i]]]$stac == 2) { @@ -503,34 +455,36 @@ DAISIE_loglik_IW <- function( if (length(status) > 0) { decstatus <- bin2dec(rev(status)) } else { - decstatus = 0 + decstatus <- 0 } - #print(status) - numcol = length(datalist) - 1 - loglik = loglik - (lgamma(M + 1) - lgamma(M - numcol + 1) - lgamma(nonendemic2 + 1) - lgamma(endemic + 1)) - #print(loglik + log(probs[1,,1,])) - loglik = loglik + log(probs[1 + nonendemic1,1 + nonendemic2,1 + endemic,1 + decstatus]) + loglik <- loglik + log(probs[1,1,1 + decstatus]) if(cond > 0) { - sysdim = 1 - totdim = lxm1 * lxm2 * lxe * sysdim - dime = list(lxm1 = lxm1,lxm2 = lxm2,lxe = lxe,sysdim = sysdim) - probs = rep(0,totdim) - probs[1] = 1 - kmi = list(kmin = 0,ki = NULL) - nndd = nndivdep(lxm2,lxe,sysdim,Kprime,M,k = 0) - parslist = list(pars = pars1,k = k,ddep = ddep,dime = dime,kmi = kmi,nndd = nndd) + sysdim <- 1 + totdim <- lxm * lxe * sysdim + dime <- list(lxm = lxm,lxe = lxe,sysdim = sysdim) + probs <- rep(0,totdim) + probs[1] <- 1 + l0ki <- list(l0 = 0,ki = NULL) + nndd <- nndivdep(lxm = lxm,lxe = lxe,sysdim = sysdim,Kprime = Kprime,M = M,k = 0,l0 = l0ki$l0) + parslist <- list(pars = pars1,k = 0,ddep = ddep,dime = dime,l0ki = l0ki,nndd = nndd) + iw_parms <- DAISIE_IW_pars(parslist) if (startsWith(methode, "odeint::")) { - probs = .Call("daisie_odeint_iw", probs, brts[(k + 1):(k + 2)], DAISIE_odeint_iw_pars(parslist), methode, abstolint, reltolint) - } - else { - y = deSolve::ode(y = probs,times = brts[(k + 1):(k + 2)],func = DAISIE_loglik_rhs_IW,parms = parslist,rtol = reltolint,atol = abstolint,method = methode) - probs = y[2,2:(totdim + 1)] + probs <- .Call("daisie_odeint_iw", probs, c(min(brts),0), iw_parms, methode, abstolint, reltolint) + } else { + y <- deSolve::ode(y = probs, + times = c(min(brts),0), + func = DAISIE_loglik_rhs_IW, + parms = iw_parms, + rtol = reltolint, + atol = abstolint, + method = methode) + probs <- y[2,2:(totdim + 1)] } - dim(probs) = c(lxm1,lxm2,lxe,sysdim) - logcond = log(1 - probs[1,1,1,1]) - loglik = loglik - logcond + dim(probs) <- c(lxm, lxe, sysdim) + logcond <- log1p(-probs[1,1,1]) + loglik <- loglik - logcond } if (pars2[4] >= 1) { s1 <- sprintf("Parameters: %f %f %f %f %f", pars1[1], pars1[2], pars1[3], pars1[4], pars1[5]) diff --git a/R/DAISIE_loglik_IW0.R b/R/DAISIE_loglik_IW0.R index 34ea2cc6..b4bb5de2 100644 --- a/R/DAISIE_loglik_IW0.R +++ b/R/DAISIE_loglik_IW0.R @@ -1,392 +1,46 @@ -kmini0 <- function(dec2binmatk, lxm, lxe, sysdim = dim(dec2binmatk)[1]) { - posc <- Matrix::rowSums(dec2binmatk) - negc <- log2(sysdim) - posc - kmin <- rep(negc, each = lxm * lxe) - dim(kmin) <- c(lxm, lxe, sysdim) - ki <- kimat(dec2binmatk) - res <- list(kmin = kmin, ki = ki) - return(res) -} - -nndivdep0 <- function(lxm, lxe, sysdim, Kprime, M, k) { - nnm <- c(0, 0:(lxm + 1)) - nne <- c(0, 0, 0:(lxe + 1)) - lnnm <- length(nnm) - lnne <- length(nne) - nn <- matrix(nnm, nrow = lnnm, ncol = lnne, byrow = F) + - matrix(nne, nrow = lnnm, ncol = lnne, byrow = T) - nn <- replicate(sysdim, nn) - nil2lxm <- 2:(lxm + 1) - nil2lxe <- 3:(lxe + 2) - nile <- rep(1, lxe) - allc <- 1:sysdim - divdepfac <- pmax(array(0, dim = c(lxm + 3, lxe + 4, sysdim)), - 1 - (nn + k) / Kprime) - divdepfacmin1 <- pmax(array(0, dim = c(lxm + 3, lxe + 4, sysdim)), - 1 - (nn + k - 1) / Kprime) - divdepfac <- divdepfac[nil2lxm, nil2lxe, allc] - divdepfacmin1 <- divdepfacmin1[nil2lxm, nil2lxe, allc] - Mminm <- M - nn[nil2lxm, nile, allc] - res <- list(nn = nn, - divdepfac = divdepfac, - divdepfacmin1 = divdepfacmin1, - Mminm = Mminm) - return(res) -} - -DAISIE_loglik_rhs_IW0 <- function(t, x, pars) { - lac <- pars[[1]][1] - mu <- pars[[1]][2] - Kprime <- pars[[1]][3] - gam <- pars[[1]][4] - laa <- pars[[1]][5] - M <- pars[[1]][6] - kk <- pars[[2]] - ddep <- pars[[3]] - lxm <- pars[[4]]$lxm - lxe <- pars[[4]]$lxe - sysdim <- pars[[4]]$sysdim - kmin <- pars[[5]]$kmin - kplus <- kk - kmin - ki <- pars[[5]]$ki - nn <- pars[[6]]$nn - divdepfac <- pars[[6]]$divdepfac - divdepfacmin1 <- pars[[6]]$divdepfacmin1 - Mminm <- pars[[6]]$Mminm - - dim(x) <- c(lxm, lxe, sysdim) - xx <- array(0, dim = c(lxm + 2, lxe + 3, sysdim)) - xx[2:(lxm + 1), 3:(lxe + 2), 1:sysdim] <- x - nil2lxm <- 2:(lxm + 1) - nil2lxe <- 3:(lxe + 2) - allc <- 1:sysdim - nile <- rep(1, lxe) - nilm <- rep(1, lxm) - Mminm2 <- Mminm + 1 - kmin - Mminm2[Mminm2 < 0] <- 0 - Mminm[Mminm < 0] <- 0 - if (sysdim == 1) { - dim(Mminm) <- c(lxm, lxe) - } - dx <- gam * divdepfacmin1 * Mminm2 * xx[nil2lxm - 1, nil2lxe, allc] + #immigration - mu * nn[nil2lxm + 1, nile, allc] * xx[nil2lxm + 1, nil2lxe, allc] + #extinction non-endemics - mu * nn[nilm, nil2lxe + 1, allc] * xx[nil2lxm, nil2lxe + 1, allc] + #extinction endemics - lac * divdepfacmin1 * nn[nil2lxm + 1, nile, allc] * xx[nil2lxm + 1, nil2lxe - 2, allc] + #cladogenesis non-endemics - lac * divdepfacmin1 * nn[nilm, nil2lxe - 1, allc] * xx[nil2lxm, nil2lxe - 1, allc] + #cladogenesis endemics - 2 * kplus * lac * divdepfacmin1 * xx[nil2lxm, nil2lxe - 1, allc] + #cladogenesis species in tree - laa * nn[nil2lxm + 1, nile, allc] * xx[nil2lxm + 1, nil2lxe - 1, allc] + #anagenesis non-endemics - -(laa * (nn[nil2lxm, nile, allc] + kmin) + (gam * divdepfac * Mminm) + - (lac * divdepfac + mu) * (nn[nil2lxm, nil2lxe, allc] + kk)) * xx[nil2lxm, nil2lxe, allc] - if (sysdim > 1) { - dx <- dx + - laa * tensor::tensor(xx[nil2lxm, nil2lxe, allc], ki, 3, 2) + # anagenesis in colonizing lineage - 2 * lac * divdepfacmin1 * tensor::tensor(xx[nil2lxm, nil2lxe - 1, allc], ki, 3, 2) # cladogenesis in colonizing lineage - } - dim(dx) <- c(sysdim * lxm * lxe, 1) - return(list(dx)) -} - - -DAISIE_loglik_IW0 <- function( - pars1, - pars2, - datalist, - methode = "ode45", - abstolint = 1E-16, - reltolint = 1E-14, - verbose = verbose) { - - brts <- c(-abs(datalist[[1]]$brts_table[, 1]), 0) - clade <- datalist[[1]]$brts_table[, 2] - event <- datalist[[1]]$brts_table[, 3] - pars1 <- as.numeric(pars1) - - ddep <- pars2[2] - cond <- pars2[3] - if (cond > 1) { - stop('cond > 1 has not been implemented in the island-wide model.') - } - - lac <- pars1[1] - mu <- pars1[2] - Kprime <- pars1[3] - if (ddep == 0) - { - Kprime <- Inf - } - gam <- pars1[4] - laa <- pars1[5] - M <- pars1[6] - - if (min(pars1) < 0) { - message("One or more parameters are negative.\n") - loglik <- -Inf - return(loglik) - } - if ((ddep == 1 | ddep == 11) & ceiling(Kprime) < length(brts)) { - if (verbose) { - cat("The proposed value of K is incompatible with the number of species - in the clade. Likelihood for this parameter set - will be set to -Inf. \n") - } - loglik <- -Inf - return(loglik) - } - if (ddep == 1 | ddep == 11) { - lx <- min(1 + ceiling(Kprime), DDD::roundn(pars2[1]) ) - } else { - lx <- DDD::roundn(pars2[1]) - } - lxm <- min(lx, M + 1) - lxe <- lx - sysdimchange <- 1 - sysdim <- 1 - totdim <- lxm * lxe * sysdim - probs <- rep(0, totdim) - probs[1] <- 1 - loglik <- 0 - expandvec <- NULL - for (k in 0:(length(brts) - 2)) { - if (pars2[4] == 2) { - print(paste("k = ", k, ", sysdim = ", sysdim, sep = "")) - utils::flush.console() - } - dime <- list(lxm = lxm, lxe = lxe, sysdim = sysdim) - if (sysdimchange == 1) { - if (sysdim > 1) { - dec2binmatk <- dec2binmat(log2(sysdim)) - kmi <- kmini0(dec2binmatk, lxm, lxe, sysdim) - } else if (sysdim == 1) { - kmi <- list(kmin = 0, ki = NULL) - } - sysdimchange <- 0 - } - nndd <- nndivdep0(lxm, lxe, sysdim, Kprime, M, k) - parslist <- list(pars = pars1, kk = k, ddep = ddep, dime = dime, kmi = kmi, nndd = nndd) - y <- deSolve::ode(y = probs, times = brts[(k + 1):(k + 2)], func = DAISIE_loglik_rhs_IW, parms = parslist, rtol = reltolint, atol = abstolint, method = methode) - probs <- y[2, 2:(totdim + 1)] - cp <- checkprobs2(NA, loglik, probs, verbose); loglik = cp[[1]]; probs = cp[[2]] - dim(probs) <- c(lxm, lxe, sysdim) - - if (k < (length(brts) - 2)) { - divdepfac <- nndd$divdepfac - if (event[k + 2] == 1) { - #Mminmminl = nndd$Mminm - kmi$kmin - Mminmminl <- nndd$Mminm - clade[k + 1] - Mminmminl[Mminmminl < 0] <- 0 - probs <- gam * divdepfac * Mminmminl * probs[ , , 1:sysdim] - probs <- c(probs, rep(0, totdim)) - - sysdim <- sysdim * 2 - expandvec <- c(expandvec, clade[k + 2]) - sysdimchange <- 1 - } else { - probs <- lac * divdepfac * probs[ , , 1:sysdim] - if (event[k + 2] == 2) { - tocollapse <- which(expandvec == clade[k + 2]) - sr <- selectrows(sysdim, 2 ^ (tocollapse - 1)) - probs <- probs[ , , sr[, 1]] + probs[ , , sr[, 2]] - sysdim <- sysdim / 2 - dim(probs) <- c(lxm, lxe, sysdim) - expandvec <- expandvec[-tocollapse] - sysdimchange <- 1 - } - } - cp <- checkprobs2(NA, loglik, probs, verbose); loglik = cp[[1]]; probs = cp[[2]] - totdim <- lxm * lxe * sysdim - dim(probs) <- c(totdim, 1) - #print(head(probs,n = 5)) - } - } - dim(probs) <- c(lxm, lxe, sysdim) - expandedclades <- which(pracma::histc(clade, 1:length(clade))$cnt == 1) - status <- rep(0, lexpandedclades <- length(expandedclades)) - if (lexpandedclades > 0) { - for (i in lexpandedclades:1) { - if (datalist[[1 + expandedclades[i]]]$stac == 2) { - status[i] <- 1 - } - } - } - endemic <- 0 - nonendemic <- 0 - for (i in 2:length(datalist)) { - endemic <- endemic + (datalist[[i]]$stac == 5) - nonendemic <- nonendemic + (datalist[[i]]$stac == 1) + (datalist[[i]]$stac == 3) - } - if (length(status) > 0) { - decstatus <- bin2dec(status) - } else { - decstatus <- 0 - } - #print(loglik + log(probs)) - loglik <- loglik + log(probs[1 + nonendemic, 1 + endemic, 1 + decstatus]) - - if (cond > 0) { - sysdim <- 1 - totdim <- lxm * lxe * sysdim - dime <- list(lxm = lxm, lxe = lxe, sysdim = sysdim) - probs <- rep(0, totdim) - probs[1] <- 1 - kmi <- list(kmin = 0, ki = NULL) - nndd <- nndivdep0(lxm, lxe, sysdim, Kprime, M, k = 0) - parslist <- list(pars = pars1, kk = k, ddep = ddep, dime = dime, kmi = kmi, nndd = nndd) - y <- deSolve::ode(y = probs, times = brts[(k + 1):(k + 2)], func = DAISIE_loglik_rhs_IW0, parms = parslist, rtol = reltolint, atol = abstolint, method = methode) - probs <- y[2, 2:(totdim + 1)] - dim(probs) <- c(lxm, lxe, sysdim) - logcond <- log(1 - probs[1, 1, 1]) - loglik <- loglik - logcond - } - - if (pars2[4] > 0) { - s1 <- sprintf("Parameters: %f %f %f %f %f %d", pars1[1], pars1[2], pars1[3], pars1[4], pars1[5], pars1[6]) - s2 <- sprintf(", Loglikelihood: %f", loglik) - cat(s1, s2, "\n", sep = "") - utils::flush.console() - } - return(as.numeric(loglik)) -} - DAISIE_loglik_IW_M1 <- function( pars1, pars2, + datalist, brts, stac, missnumspec, - methode = "ode45", + methode, abstolint = 1E-16, reltolint = 1E-14, verbose ) { - - if (is.na(pars2[4])) { - pars2[4] <- 0 - } - ddep <- pars2[2] - cond <- pars2[3] - brts <- c(-abs(brts), 0) - pars1 <- as.numeric(pars1) - lac <- pars1[1] - mu <- pars1[2] - Kprime <- pars1[3] - if (ddep == 0) { - Kprime <- Inf - } - gam <- pars1[4] - laa <- pars1[5] - pars1[6] <- 1 - M <- pars1[6] - - if (min(pars1) < 0) { - message("One or more parameters are negative.\n") - loglik <- -Inf - return(loglik) - } - if (ddep == 1 | ddep == 11) { - lx <- min(1 + ceiling(Kprime), round(pars2[1]) ) - } else { - lx <- DDD::roundn(pars2[1]) - } - lxm <- min(lx, M + 1) - lxe <- lx - sysdimchange <- 1 - sysdim <- 1 - totdim <- lxm * lxe * sysdim - probs <- rep(0, totdim) - probs[1] <- 1 - loglik <- 0 - expandvec <- NULL - for (k in 0:(length(brts) - 2)) { - if (pars2[4] == 2) { - print(paste("k = ", k, ", sysdim = ", sysdim, sep = "")) - utils::flush.console() - } - dime <- list(lxm = lxm, lxe = lxe, sysdim = sysdim) - if (sysdimchange == 1) { - if (sysdim > 1) { - dec2binmatk <- dec2binmat(log2(sysdim)) - kmi <- kmini0(dec2binmatk, lxm, lxe, sysdim) - } else if (sysdim == 1) { - kmi <- list(kmin = 0, ki = NULL) - } - sysdimchange <- 0 - } - nndd <- nndivdep0(lxm, lxe, sysdim, Kprime, M, k) - parslist <- list(pars = pars1, kk = k, ddep = ddep, dime = dime, kmi = kmi, nndd = nndd) - y <- deSolve::ode(y = probs, times = brts[(k + 1):(k + 2)], func = DAISIE_loglik_rhs_IW0, parms = parslist, rtol = reltolint, atol = abstolint, method = methode) - probs <- y[2, 2:(totdim + 1)] - cp <- checkprobs2(NA, loglik, probs, verbose); loglik = cp[[1]]; probs = cp[[2]] - dim(probs) <- c(lxm, lxe, sysdim) - - if (k < (length(brts) - 2)) { - divdepfac <- nndd$divdepfac - if(k == 0) { - #Mminmminl = nndd$Mminm - kmi$kmin - Mminmminl <- nndd$Mminm - Mminmminl[Mminmminl < 0] <- 0 - probs <- gam * divdepfac * Mminmminl * probs[ , , 1:sysdim] - probs <- c(probs, rep(0, totdim)) - sysdim <- 2 - sysdimchange <- 1 - } else { - probs <- lac * divdepfac * probs[ , , 1:sysdim] - if (k == 1) { - probs <- probs[ , , 1] + probs[ , , 2] - sysdim <- 1 - dim(probs) <- c(lxm, lxe, sysdim) - sysdimchange <- 1 - } - } - cp <- checkprobs2(NA, loglik, probs, verbose); loglik = cp[[1]]; probs = cp[[2]] - totdim <- lxm * lxe * sysdim - dim(probs) <- c(totdim, 1) - } - } - dim(probs) <- c(lxm, lxe, sysdim) - endemic <- (stac == 5) - nonendemic <- (stac == 1) + (stac == 3) - decstatus <- (stac == 2) * (sysdim > 1) #when stac = 4, decstatus = 0 - #print(probs) - loglik <- loglik + log(probs[1 + nonendemic, 1 + endemic, 1 + decstatus]) - - if (pars2[4] > 0) { - s1 <- sprintf("Parameters: %f %f %f %f %f %d", pars1[1], pars1[2], pars1[3], pars1[4], pars1[5], pars1[6]) - s2 <- sprintf(", Loglikelihood: %f", loglik) - cat(s1, s2, "\n", sep = "") - utils::flush.console() - } - return(as.numeric(loglik)) -} - -DAISIE_loglik_IW0_choosepar <- function( - trparsopt, - trparsfix, - idparsopt, - idparsfix, - pars2, - datalist, - methode, - abstolint, - reltolint -) { - trpars1 <- rep(0, 6) - trpars1[idparsopt] <- trparsopt - if (length(idparsfix) != 0) { - trpars1[idparsfix] <- trparsfix - } - if (max(trpars1) > 1 | min(trpars1) < 0) { - loglik <- -Inf + if(missnumspec > 0) { + stop('This likelihood computation cannot deal with missing species.') + } + if(!(stac %in% c(0,2,4))) { + stop('This likelihood computation must have explicit colonization times or none at all.') + } + datalist2 <- list() + datalist2[[1]] <- list(island_age = max(abs(brts)), not_present = as.numeric(is.null(datalist))) + if(!is.null(datalist)) { + pars2[3] <- 0 + datalist2[[2]] <- datalist + loglik <- DAISIE_loglik_IW( + pars1 = pars1, + pars2 = pars2, + datalist = datalist2, + methode = methode, + abstolint = abstolint, + reltolint = reltolint, + verbose = verbose) } else { - pars1 <- trpars1 / (1 - trpars1) - if (min(pars1) < 0) { - loglik <- -Inf - } else { - loglik <- DAISIE_loglik_IW0(pars1, pars2, datalist, methode, abstolint, reltolint) - } - if (is.nan(loglik) || is.na(loglik)) { - cat("There are parameter values used which cause numerical problems.\n") - loglik <- -Inf - } + loglik <- DAISIE_loglik( + pars1 = pars1, + pars2 = pars2, + brts = max(abs(brts)), + stac = 0, + missnumspec = missnumspec, + methode = 'ode45', + abstolint = abstolint, + reltolint = reltolint, + verbose = verbose + ) } return(loglik) } diff --git a/R/DAISIE_make_archipelago.R b/R/DAISIE_make_archipelago.R index 87f7efd4..07189041 100644 --- a/R/DAISIE_make_archipelago.R +++ b/R/DAISIE_make_archipelago.R @@ -35,7 +35,7 @@ DAISIE_make_archipelago <- function(archipelago, archipelago_daisie[[1]]$distance_nearest_big <- distance_nearest_big archipelago_daisie[[1]]$name <- archipelago - archipelago_daisie <- Add_brt_table(archipelago_daisie) + archipelago_daisie <- add_brt_table(archipelago_daisie) archipelago_daisie[[1]]$brts_table <- NULL return(archipelago_daisie) diff --git a/R/default_params_doc.R b/R/default_params_doc.R index ff15dbd2..bc9a577d 100644 --- a/R/default_params_doc.R +++ b/R/default_params_doc.R @@ -405,8 +405,10 @@ #' parameters being optimized by the specified amount which should be very #' small, e.g. 1e-5. Jitter when \code{link[subplex]{subplex}()} produces #' incorrect output due to parameter transformation. +#' @param num_cycles The number of cycles the optimizer will go through. +#' Default is 1. #' @param trait_pars A named list containing diversification rates considering -#' two trait states created by \code{\link{create_trait_pars}}: +#' 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} @@ -431,10 +433,10 @@ #' \itemize{ #' \item{model: the CS model to run, options are \code{1} for single rate #' DAISIE model, \code{2} for multi-rate DAISIE, or \code{0} for IW test -#' model} +#' model.} #' \item{relaxed_par: the parameter to relax (integrate over). Options are #' \code{"cladogenesis"}, \code{"extinction"}, \code{"carrying_capacity"}, -#' \code{"immigration"}, or \code{"anagenesis"}} +#' \code{"immigration"}, or \code{"anagenesis"}.} #' } #' @param DAISIE_par A numeric parameter to evaluate the integral of the #' function. @@ -561,6 +563,7 @@ default_params_doc <- function( proptime_max, current_area, jitter, + num_cycles, trait_pars, relaxed_par, relaxed_rate_pars, diff --git a/README.md b/README.md index ebf07b7d..f8d09656 100644 --- a/README.md +++ b/README.md @@ -11,6 +11,15 @@ Branch|![GHA logo](pics/github_actions_logo.png)|[![Codecov logo](pics/Codecov.p `master`|[![R build status](https://github.com/rsetienne/DAISIE/workflows/R-CMD-check/badge.svg?branch=master)](https://github.com/rsetienne/DAISIE/actions)|[![codecov.io](https://codecov.io/github/rsetienne/DAISIE/coverage.svg?branch=master)](https://codecov.io/github/rsetienne/DAISIE/branch/master) `develop`|[![R build status](https://github.com/rsetienne/DAISIE/workflows/R-CMD-check/badge.svg?branch=develop)](https://github.com/rsetienne/DAISIE/actions)|[![codecov.io](https://codecov.io/github/rsetienne/DAISIE/coverage.svg?branch=develop)](https://codecov.io/github/rsetienne/DAISIE/branch/develop) + +DAISIE is an `R` package that simulates and computes the (maximum) likelihood of a dynamical model of island biota assembly through speciation, immigration and extinction. + +The model can be fitted to both empirical dated phylogenies and simulated data. + +* For an overview of the simulation functionality see [here](https://cran.r-project.org/web/packages/DAISIE/vignettes/demo_sim.html). +* Details and an overview of the maximum likelihood inference capabilities to estimate parameters see [here](https://cran.r-project.org/web/packages/DAISIE/vignettes/demo_optimize.html). +* For details on comparing between two diversity depedence models see [here](https://cran.r-project.org/web/packages/DAISIE/vignettes/demo_CSvsIW.html). + ## Installing DAISIE **N.B.: MacOS users may experience issues when installing DAISIE, especially when on MacOS Big Sur. If that is you case, please see [here](doc/DAISIE_macOS.md) for detailed installation instructions.** @@ -73,7 +82,7 @@ Remotes: ## 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 3.2.0. 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 @@ -90,3 +99,5 @@ Valente, L., Etienne, R.S., & Garcia-R., J.C. (2019). Deep Macroevolutionary Imp 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 + +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 diff --git a/man/DAISIE_ML.Rd b/man/DAISIE_ML.Rd index 24e2315b..ac453bb7 100644 --- a/man/DAISIE_ML.Rd +++ b/man/DAISIE_ML.Rd @@ -29,7 +29,8 @@ DAISIE_ML_CS( CS_version = 1, verbose = 0, tolint = c(1e-16, 1e-10), - jitter = 0 + jitter = 0, + num_cycles = 1 ) } \arguments{ @@ -164,10 +165,10 @@ model, for a relaxed-rate model a list with the following elements: \itemize{ \item{model: the CS model to run, options are \code{1} for single rate DAISIE model, \code{2} for multi-rate DAISIE, or \code{0} for IW test - model} + model.} \item{relaxed_par: the parameter to relax (integrate over). Options are \code{"cladogenesis"}, \code{"extinction"}, \code{"carrying_capacity"}, -\code{"immigration"}, or \code{"anagenesis"}} +\code{"immigration"}, or \code{"anagenesis"}.} }} \item{verbose}{In simulation and dataprep functions a logical, @@ -184,6 +185,9 @@ tolerance of the integration.} parameters being optimized by the specified amount which should be very small, e.g. 1e-5. Jitter when \code{link[subplex]{subplex}()} produces incorrect output due to parameter transformation.} + +\item{num_cycles}{The number of cycles the optimizer will go through. +Default is 1.} } \value{ The output is a dataframe containing estimated parameters and diff --git a/man/DAISIE_ML1.Rd b/man/DAISIE_ML1.Rd index 52b6cf82..f9bf7c81 100644 --- a/man/DAISIE_ML1.Rd +++ b/man/DAISIE_ML1.Rd @@ -25,7 +25,8 @@ DAISIE_ML1( verbose = 0, tolint = c(1e-16, 1e-10), island_ontogeny = NA, - jitter = 0 + jitter = 0, + num_cycles = 1 ) } \arguments{ @@ -138,10 +139,10 @@ model, for a relaxed-rate model a list with the following elements: \itemize{ \item{model: the CS model to run, options are \code{1} for single rate DAISIE model, \code{2} for multi-rate DAISIE, or \code{0} for IW test - model} + model.} \item{relaxed_par: the parameter to relax (integrate over). Options are \code{"cladogenesis"}, \code{"extinction"}, \code{"carrying_capacity"}, -\code{"immigration"}, or \code{"anagenesis"}} +\code{"immigration"}, or \code{"anagenesis"}.} }} \item{verbose}{In simulation and dataprep functions a logical, @@ -167,6 +168,9 @@ functions \code{island_ontogeny = NA} assumes constant ontogeny.} parameters being optimized by the specified amount which should be very small, e.g. 1e-5. Jitter when \code{link[subplex]{subplex}()} produces incorrect output due to parameter transformation.} + +\item{num_cycles}{The number of cycles the optimizer will go through. +Default is 1.} } \value{ The output is a dataframe containing estimated parameters and diff --git a/man/DAISIE_ML2.Rd b/man/DAISIE_ML2.Rd index 24053218..bb54557a 100644 --- a/man/DAISIE_ML2.Rd +++ b/man/DAISIE_ML2.Rd @@ -21,7 +21,8 @@ DAISIE_ML2( optimmethod = "subplex", verbose = 0, tolint = c(1e-16, 1e-10), - jitter = 0 + jitter = 0, + num_cycles = 1 ) } \arguments{ @@ -138,6 +139,9 @@ tolerance of the integration.} parameters being optimized by the specified amount which should be very small, e.g. 1e-5. Jitter when \code{link[subplex]{subplex}()} produces incorrect output due to parameter transformation.} + +\item{num_cycles}{The number of cycles the optimizer will go through. +Default is 1.} } \value{ The output is a dataframe containing estimated parameters and diff --git a/man/DAISIE_ML3.Rd b/man/DAISIE_ML3.Rd index 2a78bd89..f39befd1 100644 --- a/man/DAISIE_ML3.Rd +++ b/man/DAISIE_ML3.Rd @@ -22,7 +22,8 @@ DAISIE_ML3( CS_version = 1, verbose = 0, tolint = c(1e-16, 1e-10), - jitter = 0 + jitter = 0, + num_cycles = 1 ) } \arguments{ @@ -121,10 +122,10 @@ model, for a relaxed-rate model a list with the following elements: \itemize{ \item{model: the CS model to run, options are \code{1} for single rate DAISIE model, \code{2} for multi-rate DAISIE, or \code{0} for IW test - model} + model.} \item{relaxed_par: the parameter to relax (integrate over). Options are \code{"cladogenesis"}, \code{"extinction"}, \code{"carrying_capacity"}, -\code{"immigration"}, or \code{"anagenesis"}} +\code{"immigration"}, or \code{"anagenesis"}.} }} \item{verbose}{In simulation and dataprep functions a logical, @@ -141,6 +142,9 @@ tolerance of the integration.} parameters being optimized by the specified amount which should be very small, e.g. 1e-5. Jitter when \code{link[subplex]{subplex}()} produces incorrect output due to parameter transformation.} + +\item{num_cycles}{The number of cycles the optimizer will go through. +Default is 1.} } \value{ The output is a dataframe containing estimated parameters and diff --git a/man/DAISIE_ML4.Rd b/man/DAISIE_ML4.Rd index f24e6c44..12ccece6 100644 --- a/man/DAISIE_ML4.Rd +++ b/man/DAISIE_ML4.Rd @@ -22,7 +22,8 @@ DAISIE_ML4( verbose = 0, tolint = c(1e-16, 1e-10), island_ontogeny = NA, - jitter = 0 + jitter = 0, + num_cycles = 1 ) } \arguments{ @@ -112,10 +113,10 @@ model, for a relaxed-rate model a list with the following elements: \itemize{ \item{model: the CS model to run, options are \code{1} for single rate DAISIE model, \code{2} for multi-rate DAISIE, or \code{0} for IW test - model} + model.} \item{relaxed_par: the parameter to relax (integrate over). Options are \code{"cladogenesis"}, \code{"extinction"}, \code{"carrying_capacity"}, -\code{"immigration"}, or \code{"anagenesis"}} +\code{"immigration"}, or \code{"anagenesis"}.} }} \item{verbose}{In simulation and dataprep functions a logical, @@ -141,6 +142,9 @@ functions \code{island_ontogeny = NA} assumes constant ontogeny.} parameters being optimized by the specified amount which should be very small, e.g. 1e-5. Jitter when \code{link[subplex]{subplex}()} produces incorrect output due to parameter transformation.} + +\item{num_cycles}{The number of cycles the optimizer will go through. +Default is 1.} } \value{ The output is a dataframe containing estimated parameters and @@ -160,7 +164,7 @@ deviation for the parameter which is allowed to vary} \item{loglik}{ gives the maximum loglikelihood} \item{df}{ gives the number -of estimated parameters, i.e. degrees of feedom} +of estimated parameters, i.e. degrees of freedom} \item{conv}{ gives a message on convergence of optimization; conv = 0 means convergence} } diff --git a/man/DAISIE_ML_IW.Rd b/man/DAISIE_ML_IW.Rd index 977144ac..324f4437 100644 --- a/man/DAISIE_ML_IW.Rd +++ b/man/DAISIE_ML_IW.Rd @@ -20,7 +20,8 @@ DAISIE_ML_IW( optimmethod = "subplex", verbose = 0, tolint = c(1e-16, 1e-14), - jitter = 0 + jitter = 0, + num_cycles = 1 ) } \arguments{ @@ -119,6 +120,9 @@ tolerance of the integration.} parameters being optimized by the specified amount which should be very small, e.g. 1e-5. Jitter when \code{link[subplex]{subplex}()} produces incorrect output due to parameter transformation.} + +\item{num_cycles}{The number of cycles the optimizer will go through. +Default is 1.} } \value{ The output is a dataframe containing estimated parameters and diff --git a/man/DAISIE_MW_ML.Rd b/man/DAISIE_MW_ML.Rd index 527d0262..3ac9ad88 100644 --- a/man/DAISIE_MW_ML.Rd +++ b/man/DAISIE_MW_ML.Rd @@ -26,7 +26,8 @@ DAISIE_MW_ML( distance_type = "continent", distance_dep = "power", parallel = "local", - cpus = 3 + cpus = 3, + num_cycles = 1 ) } \arguments{ @@ -92,8 +93,16 @@ Default is "lsodes"} "subplex" (see subplex package). Alternative is 'simplex' which was the method in previous versions.} -\item{CS_version}{For internal testing purposes only. Default is 1, the -original DAISIE code.} +\item{CS_version}{a numeric or list. Default is 1 for the standard DAISIE +model, for a relaxed-rate model a list with the following elements: +\itemize{ + \item{model: the CS model to run, options are \code{1} for single rate + DAISIE model, \code{2} for multi-rate DAISIE, or \code{0} for IW test + 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{verbose}{sets whether parameters and likelihood should be printed (1) or not (0)} @@ -135,6 +144,9 @@ machine.} \item{cpus}{Number of cpus used in parallel computing. Default is 3. Will not have an effect if parallel = 'no'.} + +\item{num_cycles}{The number of cycles the optimizer will go through. +Default is 1.} } \value{ The output is a dataframe containing estimated parameters and diff --git a/man/DAISIE_SR_ML.Rd b/man/DAISIE_SR_ML.Rd index 0deb5f18..6ce67053 100644 --- a/man/DAISIE_SR_ML.Rd +++ b/man/DAISIE_SR_ML.Rd @@ -24,7 +24,8 @@ DAISIE_SR_ML_CS( CS_version = 1, verbose = 0, tolint = c(1e-16, 1e-10), - jitter = 0 + jitter = 0, + num_cycles = 1 ) } \arguments{ @@ -70,7 +71,7 @@ time of shift} e.g. c(1,3) if lambda^c and K should not be optimized.} \item{idparsnoshift}{The ids of the parameters that should not be different -before and after the shift} +before and after the shift.} \item{res}{Sets the maximum number of species for which a probability must be computed, must be larger than the size of the largest clade} @@ -80,19 +81,20 @@ no diversity dependence \cr ddmodel = 1 : linear dependence in speciation rate \cr ddmodel = 11: linear dependence in speciation rate and in immigration rate \cr ddmodel = 2 : exponential dependence in speciation rate\cr ddmodel = 21: exponential dependence in speciation rate and in -immigration rate\cr} +immigration rate.\cr} \item{cond}{cond = 0 : conditioning on island age \cr cond = 1 : conditioning on island age and non-extinction of the island biota \cr} -\item{island_ontogeny}{type of island ontonogeny. If NA, then constant ontogeny is assumed} +\item{island_ontogeny}{type of island ontonogeny. If NA, then constant +ontogeny is assumed.} \item{tol}{Sets the tolerances in the optimization. Consists of: \cr reltolx = relative tolerance of parameter values in optimization \cr reltolf = relative tolerance of function value in optimization \cr abstolx = absolute -tolerance of parameter values in optimization} +tolerance of parameter values in optimization.} -\item{maxiter}{Sets the maximum number of iterations in the optimization} +\item{maxiter}{Sets the maximum number of iterations in the optimization.} \item{methode}{Method of the ODE-solver. See package deSolve for details. Default is "lsodes"} @@ -101,19 +103,30 @@ Default is "lsodes"} "subplex" (see subplex package). Alternative is 'simplex' which was the method in previous versions.} -\item{CS_version}{For internal testing purposes only. Default is 1, the -original DAISIE code.} +\item{CS_version}{a numeric or list. Default is 1 for the standard DAISIE +model, for a relaxed-rate model a list with the following elements: +\itemize{ + \item{model: the CS model to run, options are \code{1} for single rate + DAISIE model, \code{2} for multi-rate DAISIE, or \code{0} for IW test + 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{verbose}{sets whether parameters and likelihood should be printed (1) -or not (0)} +or not (0).} \item{tolint}{Vector of two elements containing the absolute and relative -tolerance of the integration} +tolerance of the integration.} \item{jitter}{Numeric for \code{\link[DDD]{optimizer}()}. Jitters the parameters being optimized by the specified amount which should be very small, e.g. 1e-5. Jitter when \code{link[subplex]{subplex}()} produces incorrect output due to parameter transformation.} + +\item{num_cycles}{The number of cycles the optimizer will go through. +Default is 1.} } \value{ The output is a dataframe containing estimated parameters and diff --git a/man/DAISIE_SR_loglik_CS.Rd b/man/DAISIE_SR_loglik_CS.Rd index e14d06d7..0aeb5eb7 100644 --- a/man/DAISIE_SR_loglik_CS.Rd +++ b/man/DAISIE_SR_loglik_CS.Rd @@ -90,8 +90,16 @@ applicable for endemic clades) \cr} \item{methode}{Method of the ODE-solver. See package deSolve for details. Default is "lsodes"} -\item{CS_version}{For internal testing purposes only. Default is 1, the -original DAISIE code.} +\item{CS_version}{a numeric or list. Default is 1 for the standard DAISIE +model, for a relaxed-rate model a list with the following elements: +\itemize{ + \item{model: the CS model to run, options are \code{1} for single rate + DAISIE model, \code{2} for multi-rate DAISIE, or \code{0} for IW test + 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{abstolint}{Absolute tolerance of the integration} diff --git a/man/DAISIE_convert_to_classic_plot.Rd b/man/DAISIE_convert_to_classic_plot.Rd index 702a43e1..01a2872d 100644 --- a/man/DAISIE_convert_to_classic_plot.Rd +++ b/man/DAISIE_convert_to_classic_plot.Rd @@ -11,7 +11,7 @@ DAISIE_convert_to_classic_plot(simulation_outputs, trait_pars = NULL) produced by DAISIE_sim functions.} \item{trait_pars}{A named list containing diversification rates considering -two trait states created by \code{\link{create_trait_pars}}: + 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} diff --git a/man/DAISIE_create_island.Rd b/man/DAISIE_create_island.Rd index d717259c..2bc9cbe8 100644 --- a/man/DAISIE_create_island.Rd +++ b/man/DAISIE_create_island.Rd @@ -27,7 +27,7 @@ If using an island-wide diversity dependence, this value is set to the number of mainland species.} \item{trait_pars}{A named list containing diversification rates considering -two trait states created by \code{\link{create_trait_pars}}: + 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} diff --git a/man/DAISIE_extract_stt_median.Rd b/man/DAISIE_extract_stt_median.Rd index db771784..ca08f8c5 100644 --- a/man/DAISIE_extract_stt_median.Rd +++ b/man/DAISIE_extract_stt_median.Rd @@ -18,7 +18,7 @@ with the elements \code{island_age}, \code{not_present} and \code{stt_all}. \code{nI}, \code{nA}, \code{nC} and \code{present}.} \item{trait_pars}{A named list containing diversification rates considering -two trait states created by \code{\link{create_trait_pars}}: + 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} diff --git a/man/DAISIE_format_CS.Rd b/man/DAISIE_format_CS.Rd index af7cab2e..4ee7b70c 100644 --- a/man/DAISIE_format_CS.Rd +++ b/man/DAISIE_format_CS.Rd @@ -45,7 +45,7 @@ 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}}: + 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} diff --git a/man/DAISIE_format_CS_full_stt.Rd b/man/DAISIE_format_CS_full_stt.Rd index 044bb9be..9d84dc6f 100644 --- a/man/DAISIE_format_CS_full_stt.Rd +++ b/man/DAISIE_format_CS_full_stt.Rd @@ -40,7 +40,7 @@ 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}}: + 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} diff --git a/man/DAISIE_format_CS_sampled_stt.Rd b/man/DAISIE_format_CS_sampled_stt.Rd index ec7c59a5..8920297e 100644 --- a/man/DAISIE_format_CS_sampled_stt.Rd +++ b/man/DAISIE_format_CS_sampled_stt.Rd @@ -45,7 +45,7 @@ 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}}: + 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} diff --git a/man/DAISIE_format_IW.Rd b/man/DAISIE_format_IW.Rd index 82d5f8fe..334aebe6 100644 --- a/man/DAISIE_format_IW.Rd +++ b/man/DAISIE_format_IW.Rd @@ -45,7 +45,7 @@ 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}}: + 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} diff --git a/man/DAISIE_loglik_CS.Rd b/man/DAISIE_loglik_CS.Rd index acb46a1e..a41bcb2f 100644 --- a/man/DAISIE_loglik_CS.Rd +++ b/man/DAISIE_loglik_CS.Rd @@ -84,6 +84,8 @@ the stem age of the radiation.\cr * Endemic_Singleton_MaxAge: 5 \cr * Endemic_Clade_MaxAge: 6 \cr * Endemic&Non_Endemic_Clade_MaxAge: 7 \cr \cr +* Non_endemic_MaxAge_MinAge: 8 \cr +* Endemic_Singleton_MaxAge_MinAge: 9 \cr \code{$missing_species} - number of island species that were not sampled for particular clade (only applicable for endemic clades) \cr \code{$type1or2} - whether the colonist belongs to type 1 or type 2 \cr} @@ -91,8 +93,16 @@ particular clade (only applicable for endemic clades) \cr \item{methode}{Method of the ODE-solver. See package deSolve for details. Default is "lsodes"} -\item{CS_version}{For internal testing purposes only. Default is 1, the -original DAISIE code.} +\item{CS_version}{a numeric or list. Default is 1 for the standard DAISIE +model, for a relaxed-rate model a list with the following elements: +\itemize{ + \item{model: the CS model to run, options are \code{1} for single rate + DAISIE model, \code{2} for multi-rate DAISIE, or \code{0} for IW test + 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{abstolint}{Absolute tolerance of the integration} diff --git a/man/DAISIE_loglik_integrate.Rd b/man/DAISIE_loglik_integrate.Rd index bbca7a21..3738f078 100644 --- a/man/DAISIE_loglik_integrate.Rd +++ b/man/DAISIE_loglik_integrate.Rd @@ -44,10 +44,10 @@ model, for a relaxed-rate model a list with the following elements: \itemize{ \item{model: the CS model to run, options are \code{1} for single rate DAISIE model, \code{2} for multi-rate DAISIE, or \code{0} for IW test - model} + model.} \item{relaxed_par: the parameter to relax (integrate over). Options are \code{"cladogenesis"}, \code{"extinction"}, \code{"carrying_capacity"}, -\code{"immigration"}, or \code{"anagenesis"}} +\code{"immigration"}, or \code{"anagenesis"}.} }} \item{methode}{Method of the ODE-solver. See package deSolve for details. diff --git a/man/DAISIE_plot_sims.Rd b/man/DAISIE_plot_sims.Rd index c0c22382..da1a44d9 100644 --- a/man/DAISIE_plot_sims.Rd +++ b/man/DAISIE_plot_sims.Rd @@ -36,7 +36,7 @@ divided by for plotting purposes. Larger values will lead to plots with higher resolution, but will also run slower.} \item{trait_pars}{A named list containing diversification rates considering -two trait states created by \code{\link{create_trait_pars}}: + 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} diff --git a/man/DAISIE_sim_core_trait_dependent.Rd b/man/DAISIE_sim_core_trait_dependent.Rd index 6a3fcdda..702cbcac 100644 --- a/man/DAISIE_sim_core_trait_dependent.Rd +++ b/man/DAISIE_sim_core_trait_dependent.Rd @@ -95,7 +95,7 @@ created by \code{\link{create_area_pars}()}: rate preventing it from being too large and slowing down simulation.} \item{trait_pars}{A named list containing diversification rates considering -two trait states created by \code{\link{create_trait_pars}}: + 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} diff --git a/man/DAISIE_sim_trait_dependent.Rd b/man/DAISIE_sim_trait_dependent.Rd index fc553bfc..73e81189 100644 --- a/man/DAISIE_sim_trait_dependent.Rd +++ b/man/DAISIE_sim_trait_dependent.Rd @@ -139,7 +139,7 @@ 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}}: + 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} diff --git a/man/DAISIE_sim_update_state_trait_dependent.Rd b/man/DAISIE_sim_update_state_trait_dependent.Rd index 5e01ab70..e326f285 100644 --- a/man/DAISIE_sim_update_state_trait_dependent.Rd +++ b/man/DAISIE_sim_update_state_trait_dependent.Rd @@ -32,7 +32,7 @@ of species.} \item{stt_table}{Matrix with number of species at each time step.} \item{trait_pars}{A named list containing diversification rates considering -two trait states created by \code{\link{create_trait_pars}}: + 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} diff --git a/man/are_trait_pars.Rd b/man/are_trait_pars.Rd index d75f70c8..c6e15c44 100644 --- a/man/are_trait_pars.Rd +++ b/man/are_trait_pars.Rd @@ -8,7 +8,7 @@ are_trait_pars(trait_pars) } \arguments{ \item{trait_pars}{A named list containing diversification rates considering -two trait states created by \code{\link{create_trait_pars}}: + 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} diff --git a/man/default_params_doc.Rd b/man/default_params_doc.Rd index 5a606e69..c1030959 100644 --- a/man/default_params_doc.Rd +++ b/man/default_params_doc.Rd @@ -113,6 +113,7 @@ default_params_doc( proptime_max, current_area, jitter, + num_cycles, trait_pars, relaxed_par, relaxed_rate_pars, @@ -460,10 +461,10 @@ model, for a relaxed-rate model a list with the following elements: \itemize{ \item{model: the CS model to run, options are \code{1} for single rate DAISIE model, \code{2} for multi-rate DAISIE, or \code{0} for IW test - model} + model.} \item{relaxed_par: the parameter to relax (integrate over). Options are \code{"cladogenesis"}, \code{"extinction"}, \code{"carrying_capacity"}, -\code{"immigration"}, or \code{"anagenesis"}} +\code{"immigration"}, or \code{"anagenesis"}.} }} \item{tolint}{Vector of two elements containing the absolute and relative @@ -652,8 +653,11 @@ parameters being optimized by the specified amount which should be very small, e.g. 1e-5. Jitter when \code{link[subplex]{subplex}()} produces incorrect output due to parameter transformation.} +\item{num_cycles}{The number of cycles the optimizer will go through. +Default is 1.} + \item{trait_pars}{A named list containing diversification rates considering -two trait states created by \code{\link{create_trait_pars}}: + 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} diff --git a/man/get_ana_rate.Rd b/man/get_ana_rate.Rd index af310df4..6833c95e 100644 --- a/man/get_ana_rate.Rd +++ b/man/get_ana_rate.Rd @@ -16,7 +16,7 @@ species (a.k.a non-endemic species).} of species.} \item{trait_pars}{A named list containing diversification rates considering -two trait states created by \code{\link{create_trait_pars}}: + 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} diff --git a/man/get_clado_rate.Rd b/man/get_clado_rate.Rd index 77bcf7ad..6d9fa18b 100644 --- a/man/get_clado_rate.Rd +++ b/man/get_clado_rate.Rd @@ -32,7 +32,7 @@ calculations as returned by \code{\link{create_hyper_pars}()}: \item{A}{A numeric value for island area at a given point in time.} \item{trait_pars}{A named list containing diversification rates considering -two trait states created by \code{\link{create_trait_pars}}: + 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} diff --git a/man/get_ext_rate.Rd b/man/get_ext_rate.Rd index 6efbc5f2..79a5c6de 100644 --- a/man/get_ext_rate.Rd +++ b/man/get_ext_rate.Rd @@ -33,7 +33,7 @@ rate preventing it from being too large and slowing down simulation.} \item{A}{A numeric value for island area at a given point in time.} \item{trait_pars}{A named list containing diversification rates considering -two trait states created by \code{\link{create_trait_pars}}: + 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} diff --git a/man/get_immig_rate.Rd b/man/get_immig_rate.Rd index 9f251618..d367f752 100644 --- a/man/get_immig_rate.Rd +++ b/man/get_immig_rate.Rd @@ -30,7 +30,7 @@ If using an island-wide diversity dependence, this value is set to the number of mainland species.} \item{trait_pars}{A named list containing diversification rates considering -two trait states created by \code{\link{create_trait_pars}}: + 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} diff --git a/man/update_rates.Rd b/man/update_rates.Rd index 97e28020..92cd3834 100644 --- a/man/update_rates.Rd +++ b/man/update_rates.Rd @@ -98,7 +98,7 @@ If using an island-wide diversity dependence, this value is set to the number of mainland species.} \item{trait_pars}{A named list containing diversification rates considering -two trait states created by \code{\link{create_trait_pars}}: + 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} diff --git a/src/DAISIE_IW.cpp b/src/DAISIE_IW.cpp index 3dc666f8..b29d2b33 100644 --- a/src/DAISIE_IW.cpp +++ b/src/DAISIE_IW.cpp @@ -14,6 +14,7 @@ #include #include #include +#include #include #include #include @@ -28,8 +29,17 @@ namespace { using index_v = EIGEN_DEFAULT_DENSE_INDEX_TYPE; - template - using index_t = DSizes; + + template + using index_t = std::array; + + + template + using index_array_t = std::array, J>; + + + template + using padding_t = std::array, I>; template @@ -45,20 +55,20 @@ namespace { template - index_t iofs(index_v i0, index_v i1, index_v i2); + index_t iofs(int i0, int i1); template <> - index_t<3> iofs<3>(index_v i0, index_v i1, index_v i2) + index_t<2> iofs<2>(int i0, int i1) { - return index_t<3>{i0, i1, i2}; + return index_t<2>{i0, i1}; } template <> - index_t<4> iofs<4>(index_v i0, index_v i1, index_v i2) + index_t<3> iofs<3>(int i0, int i1) { - return index_t<4>{i0, i1, i2, 0}; + return index_t<3>{i0, i1, 0}; } @@ -80,12 +90,6 @@ namespace { void rhs(const double* x, double* rdx, ThreadPoolDevice* dev); private: - static tmap mapt(List pars, const char* name) - { - DoubleVector v = pars[name]; - auto dim = dim_to_index(v); - return tmap(v.begin(), dim); - } static ctmap cmapt(List pars, const char* name) { DoubleVector v = pars[name]; @@ -94,30 +98,26 @@ namespace { } double laa_; - std::array c_; + std::array c_; matrix ki_; }; template - cpp_daisie_iw::cpp_daisie_iw(List pars) - : laa_(pars["laa"]), - c_{ - cmapt(pars, "c1"), - cmapt(pars, "c2"), - cmapt(pars, "c3"), - cmapt(pars, "c4"), - cmapt(pars, "c5"), - cmapt(pars, "c6"), - cmapt(pars, "c7"), - cmapt(pars, "c89"), - cmapt(pars, "c10"), - cmapt(pars, "c11"), - cmapt(pars, "c12"), - cmapt(pars, "c13") - } + cpp_daisie_iw::cpp_daisie_iw(List pars) : + laa_(pars["laa"]), + c_{ + cmapt(pars, "c1"), + cmapt(pars, "c2"), + cmapt(pars, "c3"), + cmapt(pars, "c4"), + cmapt(pars, "c5"), + cmapt(pars, "c6"), + cmapt(pars, "c7"), + cmapt(pars, "c8") + } { - if (rank == 4) { + if (rank > 2) { DoubleVector ki = pars["ki"]; auto dim = dim_to_index<2>(ki); ki_ = cmmap(ki.begin(), dim); @@ -125,65 +125,48 @@ namespace { } +# define xx_slice(x, y) xx.slice(iofs(x,y), dim_c) + + template <> - void cpp_daisie_iw<3>::rhs(const double* rx, double* rdx, ThreadPoolDevice* dev) + void cpp_daisie_iw<2>::rhs(const double* rx, double* rdx, ThreadPoolDevice* dev) { const auto dim_c = c_[0].dimensions(); tmap dx(rdx, dim_c); - std::array, rank> xxpad = { - std::make_pair(1,1), - std::make_pair(1,1), - std::make_pair(2,1) - }; ctmap x(rx, dim_c); - auto xx = x.pad(xxpad); - auto xce = xx.slice(iofs(1,1,1), dim_c); - auto x12 = xx.slice(iofs(1,1,2), dim_c); + auto xxpad = padding_t{{ {1,1}, {2,1} }}; + const auto xx = x.pad(xxpad); auto ddx = - c_[0] * xx.slice(iofs(0,1,2), dim_c) + - c_[1] * xx.slice(iofs(1,0,2), dim_c) + - c_[2] * xx.slice(iofs(2,1,2), dim_c) + - c_[3] * xx.slice(iofs(1,2,2), dim_c) + - c_[4] * xx.slice(iofs(1,1,3), dim_c) + - c_[5] * xx.slice(iofs(2,1,0), dim_c) + - c_[6] * xx.slice(iofs(1,2,0), dim_c) + - c_[7] * xce + - c_[8] * xx.slice(iofs(2,1,1), dim_c) + - c_[9] * xx.slice(iofs(1,2,1), dim_c) + - c_[10] * x12; + c_[0] * xx_slice(0,2) + + c_[1] * xx_slice(2,2) + + c_[2] * xx_slice(1,3) + + c_[3] * xx_slice(2,1) + + c_[4] * xx_slice(2,0) + + c_[5] * xx_slice(1,1) - + c_[6] * xx_slice(1,2); dx.device(*dev) = ddx; } template <> - void cpp_daisie_iw<4>::rhs(const double* rx, double* rdx, ThreadPoolDevice* dev) + void cpp_daisie_iw<3>::rhs(const double* rx, double* rdx, ThreadPoolDevice* dev) { const auto dim_c = c_[0].dimensions(); tmap dx(rdx, dim_c); - std::array, rank> xxpad = { - std::make_pair(1,1), - std::make_pair(1,1), - std::make_pair(2,1), - std::make_pair(0,0) - }; ctmap x(rx, dim_c); - auto xx = x.pad(xxpad); - auto xce = xx.slice(iofs(1,1,1), dim_c); - auto x12 = xx.slice(iofs(1,1,2), dim_c); + const auto product_dims = padding_t<1>{{ {2,1} }}; + auto xxpad = padding_t{{ {1,1}, {2,1}, {0,0} }}; + const auto xx = x.pad(xxpad); auto ddx = - c_[0] * xx.slice(iofs(0,1,2), dim_c) + - c_[1] * xx.slice(iofs(1,0,2), dim_c) + - c_[2] * xx.slice(iofs(2,1,2), dim_c) + - c_[3] * xx.slice(iofs(1,2,2), dim_c) + - c_[4] * xx.slice(iofs(1,1,3), dim_c) + - c_[5] * xx.slice(iofs(2,1,0), dim_c) + - c_[6] * xx.slice(iofs(1,2,0), dim_c) + - c_[7] * xce + - c_[8] * xx.slice(iofs(2,1,1), dim_c) + - c_[9] * xx.slice(iofs(1,2,1), dim_c) + - c_[10] * x12; - const array, 1> product_dims = { std::make_pair(3, 1) }; - dx.device(*dev) = ddx + (laa_ * x12 + c_[11] * xce).contract(ki_, product_dims); + c_[0] * xx_slice(0,2) + + c_[1] * xx_slice(2,2) + + c_[2] * xx_slice(1,3) + + c_[3] * xx_slice(2,1) + + c_[4] * xx_slice(2,0) + + c_[5] * xx_slice(1,1) - + c_[6] * xx_slice(1,2) + + (laa_ * xx_slice(1,2) + c_[7] * xx_slice(1,1)).contract(ki_, product_dims); + dx.device(*dev) = ddx; } @@ -192,28 +175,27 @@ namespace { std::unique_ptr pool; std::unique_ptr dev; + std::unique_ptr> iw2; std::unique_ptr> iw3; - std::unique_ptr> iw4; - daisie_iw_wrapper(List pars) + daisie_iw_wrapper(List pars) : + pool(new ThreadPool(std::thread::hardware_concurrency())), + dev(new ThreadPoolDevice(pool.get(), std::thread::hardware_concurrency())) { - pool.reset(new ThreadPool(std::thread::hardware_concurrency())); - dev.reset(new ThreadPoolDevice(pool.get(), std::thread::hardware_concurrency())); - int sysdim = pars["sysdim"]; if (1 == sysdim) { - iw3 = std::make_unique>(pars); + iw2 = std::make_unique>(pars); } else { - iw4 = std::make_unique>(pars); + iw3 = std::make_unique>(pars); } } // odeint interface void operator()(const std::vector& x, std::vector& dxdt, double) { - (iw3) ? iw3->rhs(x.data(), dxdt.data(), dev.get()) - : iw4->rhs(x.data(), dxdt.data(), dev.get()); + (iw2) ? iw2->rhs(x.data(), dxdt.data(), dev.get()) + : iw3->rhs(x.data(), dxdt.data(), dev.get()); } }; diff --git a/tests/testthat/test-DAISIE_ML1.R b/tests/testthat/test-DAISIE_ML1.R index eb8c0c11..6d0e1a8f 100644 --- a/tests/testthat/test-DAISIE_ML1.R +++ b/tests/testthat/test-DAISIE_ML1.R @@ -3,9 +3,6 @@ context("DAISIE_ML1") test_that("use", { skip_if(Sys.getenv("CI") == "", message = "Run only on CI") - # This is a rough MLE test, built for fast execution. A more thorough test - # can be found in the GitHub repository Neves-P/DAISIEtesting - utils::data(Galapagos_datalist) datalist <- Galapagos_datalist initparsopt <- c(2.5, 2.7, 20, 0.009, 1.01) @@ -13,25 +10,27 @@ test_that("use", { idparsopt <- 1:5 parsfix <- NULL idparsfix <- NULL - tested_MLE <- DAISIE:::DAISIE_ML1( - datalist = datalist, - initparsopt = initparsopt, - idparsopt = idparsopt, - parsfix = parsfix, - ddmodel = ddmodel, - idparsfix = idparsfix, - verbose = 0, - tol = c(0.01, 0.1, 0.001), - res = 15, - tolint = c(0.1, 0.01) - ) + invisible(capture.output( + tested_MLE <- DAISIE:::DAISIE_ML1( + datalist = datalist, + initparsopt = initparsopt, + idparsopt = idparsopt, + parsfix = parsfix, + ddmodel = ddmodel, + idparsfix = idparsfix, + verbose = 0, + tol = c(0.01, 0.1, 0.001), + res = 15, + tolint = c(0.1, 0.01) + ) + )) expected_MLE <- data.frame( - lambda_c = 4.1558928147225656, - mu = 4.9989880317542603, - K = 514.14066875326353, - gamma = 0.020398370167867434, - lambda_a = 3.7589983231693607, - loglik = -63.993218343764838, + lambda_c = 3.689104200780465, + mu = 4.31030299415995, + K = 906.6501180193454, + gamma = 0.0173458887696076, + lambda_a = 3.677789527566334, + loglik = -64.22005531779574, df = 5L, conv = 0L ) diff --git a/tests/testthat/test-DAISIE_ML2.R b/tests/testthat/test-DAISIE_ML2.R index 6c7705a0..8ad43afa 100644 --- a/tests/testthat/test-DAISIE_ML2.R +++ b/tests/testthat/test-DAISIE_ML2.R @@ -1,58 +1,56 @@ context("DAISIE_ML2") test_that("use", { - # This is a rough MLE test, built for fast execution. A more thorough test - # can be found in the GitHub repository Neves-P/DAISIEtesting - skip_if(Sys.getenv("CI") == "", message = "Run only on CI") utils::data(Macaronesia_datalist, package = "DAISIE") - tested_MLE <- DAISIE:::DAISIE_ML2( - datalist = Macaronesia_datalist, - initparsopt = c( - 1.053151832, - 0.052148979, - 0.512939011, - 0.133766934, - 0.152763179 - ), - idparsmat = rbind( - 1:5, - c(6, 2, 3, 7, 5), - 1:5,1:5 - ), - idparsopt = c(2, 4, 5, 6, 7), - parsfix = c(0, Inf), - idparsfix = c(1, 3), - tol = c(0.01, 0.1, 0.001), - res = 15, - tolint = c(0.1, 0.01) - ) - + invisible(capture.output( + tested_MLE <- DAISIE:::DAISIE_ML2( + datalist = Macaronesia_datalist, + initparsopt = c( + 1.053151832, + 0.052148979, + 0.512939011, + 0.133766934, + 0.152763179 + ), + idparsmat = rbind( + 1:5, + c(6, 2, 3, 7, 5), + 1:5,1:5 + ), + idparsopt = c(2, 4, 5, 6, 7), + parsfix = c(0, Inf), + idparsfix = c(1, 3), + tol = c(0.01, 0.1, 0.001), + res = 15, + tolint = c(0.1, 0.01) + ) + )) expected_MLE <- data.frame( lambda_c = c(0.0, - 0.3397037946009278, + 0.333180238362478, 0.0, 0.0), - mu = c(2.8044328568380257, - 2.8044328568380257, - 2.8044328568380257, - 2.8044328568380257), + mu = c(2.718309102127384, + 2.718309102127384, + 2.718309102127384, + 2.718309102127384), K = c(Inf, Inf, Inf, Inf), - gamma = c(0.14401353274629849, - 0.33209131444203649, - 0.14401353274629849, - 0.14401353274629849), - lambda_a = c(1.2641967865150285, - 1.2641967865150285, - 1.2641967865150285 , - 1.2641967865150285), - loglik = c(-413.80819700475951, - -413.80819700475951, - -413.80819700475951, - -413.80819700475951), + gamma = c(0.1389101810794861, + 0.3240600655694914, + 0.1389101810794861, + 0.1389101810794861), + lambda_a = c(1.17278055953365, + 1.17278055953365, + 1.17278055953365 , + 1.17278055953365), + loglik = c(-411.9377194984208, + -411.9377194984208, + -411.9377194984208, + -411.9377194984208), df = c(5L, 5L, 5L, 5L), conv = c(0L, 0L, 0L, 0L) ) diff --git a/tests/testthat/test-DAISIE_ML3.R b/tests/testthat/test-DAISIE_ML3.R index 9796ec57..41327aee 100644 --- a/tests/testthat/test-DAISIE_ML3.R +++ b/tests/testthat/test-DAISIE_ML3.R @@ -1,9 +1,8 @@ context("DAISIE_ML3") test_that("use", { skip_if(Sys.getenv("CI") == "", message = "Run only on CI") - skip("WIP") - # This is a rough MLE test, built for fast execution. A more thorough test - # can be found in the GitHub repository Neves-P/DAISIEtesting + # THIS FUNCTION DOESN'T WORK CORRECTLY YET! FOR NOW, WE TEST IT THROWS AN + # APPROPRIATE ERROR utils::data(Galapagos_datalist, package = "DAISIE") pars1 <- c(0.2, 0.1, 17, 0.001, 0.3) @@ -19,14 +18,16 @@ test_that("use", { gam = pars1[4], laa = pars1[5] ) - tested_MLE <- DAISIE:::DAISIE_ML3( + expect_error(tested_MLE <- DAISIE:::DAISIE_ML3( datalist = Galapagos_datalist, initparsopt = pars1_td[5:10], idparsopt = 5:10, parsfix = pars1_td[1:4], idparsfix = 1:4, island_ontogeny = 1 - ) + ), regexp = "This functionality is still under development and is not available yet.") + + # All code below refers to future reference test when function is completed idpars <- sort(c(5:10, 1:4)) expected_MLE <- data.frame( lambda_c = c(0.0, @@ -56,7 +57,7 @@ test_that("use", { df = c(5L, 5L, 5L, 5L), conv = c(0L, 0L, 0L, 0L) ) - expect_equal(tested_MLE, expected_MLE) + }) test_that("abuse", { diff --git a/tests/testthat/test-DAISIE_ML4.R b/tests/testthat/test-DAISIE_ML4.R index b48d8906..a7aabe31 100644 --- a/tests/testthat/test-DAISIE_ML4.R +++ b/tests/testthat/test-DAISIE_ML4.R @@ -15,20 +15,23 @@ test_that("DAISIE_ML4 is silent and produces correct output", { }) test_that("DAISIE_loglik_all_choosepar4 is silent and produces correct output", { - skip("Takes too long and produces DLSODES warnings") utils::data(Galapagos_datalist) - output <- DAISIE_loglik_all_choosepar4( - trparsopt = c(0.718364388965934, 0.728515721595379, 0.909090909090909, - 0.009245787662330, 0.502505659845103, 0.500000000000000), - trparsfix = numeric(0), - idparsopt = c(1, 2, 3, 4, 5, 6), - idparsfix = NULL, - pars2 = c(100, 0, 0, 0, NA), - datalist = Galapagos_datalist, - methode = "lsodes", - CS_version = list(model = 2, - relaxed_par = "cladogenesis"), - abstolint = 1e-16, - reltolint = 1e-10 - ) + skip_if(Sys.getenv("CI") == "", message = "Run only on CI") + # Throws warnings and DLSODES output + invisible(capture.output(suppressWarnings( + output <- DAISIE_loglik_all_choosepar4( + trparsopt = c(0.718364388965934, 0.728515721595379, 0.909090909090909, + 0.009245787662330, 0.502505659845103, 0.500000000000000), + trparsfix = numeric(0), + idparsopt = c(1, 2, 3, 4, 5, 6), + idparsfix = NULL, + pars2 = c(100, 0, 0, 0, NA), + datalist = Galapagos_datalist, + methode = "lsodes", + CS_version = list(model = 2, + relaxed_par = "cladogenesis"), + abstolint = 1e-16, + reltolint = 1e-10 + )))) + expect_equal(output, -77.50300643925) }) diff --git a/tests/testthat/test-DAISIE_ML_CS.R b/tests/testthat/test-DAISIE_ML_CS.R index e7a4a5cb..125d79ad 100644 --- a/tests/testthat/test-DAISIE_ML_CS.R +++ b/tests/testthat/test-DAISIE_ML_CS.R @@ -1,7 +1,7 @@ context("DAISIE_ML_CS") test_that("relaxed-rate DAISIE_ML_CS produces correct output", { - skip("test takes too long atm") + skip("Too slow to run") utils::data(Galapagos_datalist) CS_version <- create_CS_version(model = 2, relaxed_par = "cladogenesis") @@ -18,7 +18,7 @@ test_that("relaxed-rate DAISIE_ML_CS produces correct output", { }) test_that("relaxed-rate DAISIE_ML_CS produces correct output using simplex", { - skip("test takes too long atm") + skip("Too slow to run") utils::data(Galapagos_datalist) CS_version <- create_CS_version(model = 2, relaxed_par = "cladogenesis") diff --git a/tests/testthat/test-DAISIE_MW_ML.R b/tests/testthat/test-DAISIE_MW_ML.R index 9962bcab..718e3af9 100644 --- a/tests/testthat/test-DAISIE_MW_ML.R +++ b/tests/testthat/test-DAISIE_MW_ML.R @@ -13,40 +13,42 @@ test_that("DAISIE_MW_ML produces correct output", { 1.504297e-01, Inf, 0.000000e+00, - 6.725644e+01, + 67.2564367200001, 2.936351e-01, 5.909687e-02, 3.826885e-01, 2.651078e-02, - -3.653315e+03, + -3651.6624531664, 8.000000e+00, 0.000000e+00 ) - M19_computation <- DAISIE::DAISIE_MW_ML( - datalist = archipelagos41, - initparsopt = c( - 0.040073803, - 1.945656546, - 0.150429656, - 67.25643672, - 0.293635061, - 0.059096872, - 0.382688527, - 0.026510781), - idparsopt = c(1, 3, 4, 7, 8, 9, 10, 11), - parsfix = c(0, Inf, 0) , - idparsfix = c(2, 5, 6), - res = 100, - ddmodel = 0, - methode = 'lsodes', - cpus = 4, - parallel = 'no', - optimmethod = 'subplex', - tol = c(1E-1, 1E-3, 1E-5), - distance_type = 'continent', - distance_dep = 'area_interactive_clado' - ) + invisible(capture.output( + M19_computation <- DAISIE::DAISIE_MW_ML( + datalist = archipelagos41, + initparsopt = c( + 0.040073803, + 1.945656546, + 0.150429656, + 67.25643672, + 0.293635061, + 0.059096872, + 0.382688527, + 0.026510781), + idparsopt = c(1, 3, 4, 7, 8, 9, 10, 11), + parsfix = c(0, Inf, 0) , + idparsfix = c(2, 5, 6), + res = 100, + ddmodel = 0, + methode = 'lsodes', + cpus = 4, + parallel = 'no', + optimmethod = 'subplex', + tol = c(1E-1, 1E-3, 1E-5), + distance_type = 'continent', + distance_dep = 'area_interactive_clado' + ) + )) testthat::expect_equal( M19_Nature_parameters, @@ -68,40 +70,42 @@ test_that("DAISIE_MW_ML produces correct output when in parallel", { 1.504297e-01, Inf, 0.000000e+00, - 6.725644e+01, + 67.2564367200001, 2.936351e-01, 5.909687e-02, 3.826885e-01, 2.651078e-02, - -3.653315e+03, + -3651.6624531664, 8.000000e+00, 0.000000e+00 ) - M19_computation <- DAISIE::DAISIE_MW_ML( - datalist = archipelagos41, - initparsopt = c( - 0.040073803, - 1.945656546, - 0.150429656, - 67.25643672, - 0.293635061, - 0.059096872, - 0.382688527, - 0.026510781), - idparsopt = c(1, 3, 4, 7, 8, 9, 10, 11), - parsfix = c(0, Inf, 0) , - idparsfix = c(2, 5, 6), - res = 100, - ddmodel = 0, - methode = 'lsodes', - cpus = 4, - parallel = 'local', - optimmethod = 'subplex', - tol = c(1E-1, 1E-3, 1E-5), - distance_type = 'continent', - distance_dep = 'area_interactive_clado' - ) + invisible(capture.output( + M19_computation <- DAISIE::DAISIE_MW_ML( + datalist = archipelagos41, + initparsopt = c( + 0.040073803, + 1.945656546, + 0.150429656, + 67.25643672, + 0.293635061, + 0.059096872, + 0.382688527, + 0.026510781), + idparsopt = c(1, 3, 4, 7, 8, 9, 10, 11), + parsfix = c(0, Inf, 0) , + idparsfix = c(2, 5, 6), + res = 100, + ddmodel = 0, + methode = 'lsodes', + cpus = 4, + parallel = 'local', + optimmethod = 'subplex', + tol = c(1E-1, 1E-3, 1E-5), + distance_type = 'continent', + distance_dep = 'area_interactive_clado' + ) + )) testthat::expect_equal( M19_Nature_parameters, diff --git a/tests/testthat/test-DAISIE_SR_loglik.R b/tests/testthat/test-DAISIE_SR_loglik.R index 8a3f4fa4..467ebfa7 100644 --- a/tests/testthat/test-DAISIE_SR_loglik.R +++ b/tests/testthat/test-DAISIE_SR_loglik.R @@ -1,6 +1,6 @@ context("DAISIE_SR_loglik") -test_that("The SR simulation and inference code works", { +test_that("The SR loglik code works", { Biwa_datalist <- NULL rm(Biwa_datalist) utils::data(Biwa_datalist, package = "DAISIE") @@ -21,25 +21,7 @@ test_that("The SR simulation and inference code works", { loglik <- DAISIE_SR_loglik_CS(pars1 = pars1, pars2 = pars2, datalist = Biwa_datalist) - testthat::expect_equal(loglik, -232.0618, tol = 1E-3) - - # Simulate fish diversity over 4 Ma - set.seed(1) - M <- 312 - IslandAge <- 4 - sims <- DAISIE_SR_sim( - time = 4, - M = M - 17, - pars = pars1, - replicates = 1, - plot_sims = FALSE, - ddep = 11, - verbose = FALSE - ) - # Compare richnesses of the last time bin - testthat::expect_equal( - unname(sims[[1]][[1]]$stt_all[26, ]), c(0, 56, 11, 0, 66) - ) + testthat::expect_equal(loglik, -227.7108, tol = 1E-3) Macaronesia_datalist <- NULL rm(Macaronesia_datalist) @@ -60,13 +42,13 @@ test_that("The SR simulation and inference code works", { pars1mat <- matrix(pars1, nrow = 8, ncol = 11, byrow = T) expected_loglik <- c( -Inf, - -266.9125, + -257.7311, -Inf, - -262.7665, + -252.5254, -Inf, - -291.6965, + -280.0562, -Inf, - -261.9222 + -252.6193 ) # No shift in cladogenesis rate older before # colonization of diversifying lineages diff --git a/tests/testthat/test-DAISIE_format_CS.R b/tests/testthat/test-DAISIE_format_CS.R index e9ba5716..b3cf153b 100644 --- a/tests/testthat/test-DAISIE_format_CS.R +++ b/tests/testthat/test-DAISIE_format_CS.R @@ -135,46 +135,6 @@ test_that("output with empty island and verbose = TRUE", { ) }) -test_that("silent with empty 2 type island", { - skip("this scenario can't run (min_type2 always produces species)") - pars <- c(0, 10, 1, 0.0001, 0, 0, 10, 1, 0.0001, 0) - totaltime <- 1 - M <- 1 - mainland_n <- M - replicates <- 1 - prop_type2_pool <- 0.1 - verbose <- FALSE - sample_freq <- 1 - set.seed(1) - island_replicates <- list() - island_replicates <- DAISIE:::DAISIE_sim_min_type2( - time = totaltime, - M = M, - pars = pars, - replicates = replicates, - prop_type2_pool = prop_type2_pool, - verbose = verbose, - 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) - ) - expect_silent( - formatted_CS_sim <- DAISIE:::DAISIE_format_CS( # sim_min_type2 produces list with one extra element and fails - island_replicates = island_replicates, - time = totaltime, - M = mainland_n, - sample_freq = sample_freq, - verbose = verbose - ) - ) -}) - test_that("silent with non-empty 2 type island", { pars <- c(0.4, 0.1, 10, 1, 0.5, 0.4, 0.1, 10, 1, 0.5) totaltime <- 1 diff --git a/tests/testthat/test-DAISIE_format_IW.R b/tests/testthat/test-DAISIE_format_IW.R index 73676f3b..390b84d5 100644 --- a/tests/testthat/test-DAISIE_format_IW.R +++ b/tests/testthat/test-DAISIE_format_IW.R @@ -42,14 +42,14 @@ test_that("silent with empty island with correct output", { colnames(stt_all) <- c("Time", "nI", "nA", "nC") stt_all[1, ] <- c(1, 0, 0, 0) stt_all[2, ] <- c(0, 0, 0, 0) - brts_table <- matrix(ncol = 4, nrow = 1) - colnames(brts_table) <- c("brt", "clade", "event", "endemic") - brts_table[1, ] <- c(1, 0, 0, NA) + brts_table <- matrix(ncol = 5, nrow = 1) + colnames(brts_table) <- c("brt", "clade", "event", "endemic", "col") + brts_table[1, ] <- c(1, 0, 0, NA, NA) expected_IW_format[[1]][[1]] <- list(island_age = 1, not_present = 10, stt_all = stt_all, brts_table = brts_table) - expect_true(all.equal(formated_IW_sim, expected_IW_format, tolerance = 1e-7)) + expect_equal(formated_IW_sim, expected_IW_format, tolerance = 1e-7) }) test_that("silent with non-empty island with correct output", { @@ -93,14 +93,14 @@ test_that("silent with non-empty island with correct output", { colnames(stt_all) <- c("Time", "nI", "nA", "nC") stt_all[1, ] <- c(1, 0, 0, 0) stt_all[2, ] <- c(0, 2, 0, 3) - brts_table <- matrix(ncol = 4, nrow = 6) - colnames(brts_table) <- c("brt", "clade", "event", "endemic") - brts_table[1, ] <- c(1, 0, 0, NA) - brts_table[2, ] <- c(0.9244818166871660, 1, 1, 1) - brts_table[3, ] <- c(0.9105856673960619, 1, 2, 1) - brts_table[4, ] <- c(0.5557734125062590, 2, 1, 0) - brts_table[5, ] <- c(0.5288428248966160, 3, 1, 0) - brts_table[6, ] <- c(0.3146835586399670, 1, 3, 1) + brts_table <- matrix(ncol = 5, nrow = 6) + colnames(brts_table) <- c("brt", "clade", "event", "endemic", "col") + brts_table[1, ] <- c(1, 0, 0, NA, NA) + brts_table[2, ] <- c(0.9244818166871660, 1, 1, 1, NA) + brts_table[3, ] <- c(0.9105856673960619, 1, 2, 1, NA) + brts_table[4, ] <- c(0.5557734125062590, 2, 1, 0, NA) + brts_table[5, ] <- c(0.5288428248966160, 3, 1, 0, NA) + brts_table[6, ] <- c(0.3146835586399670, 1, 3, 1, NA) expected_IW_format[[1]][[1]] <- list(island_age = 1, not_present = 7, stt_all = stt_all, @@ -119,7 +119,7 @@ test_that("silent with non-empty island with correct output", { 0.5288428248966160), stac = 4, missing_species = 0) - expect_true(all.equal(formated_IW_sim, expected_IW_format, tolerance = 1e-7)) + expect_equal(formated_IW_sim, expected_IW_format, tolerance = 1e-7) }) test_that("DAISIE_format_IW prints when verbose = TRUE", { @@ -201,9 +201,9 @@ test_that("silent with empty nonoceanic island with correct output", { colnames(stt_all) <- c("Time", "nI", "nA", "nC") stt_all[1, ] <- c(1, 1, 2, 0) stt_all[2, ] <- c(0, 0, 0, 0) - brts_table <- matrix(ncol = 4, nrow = 1) - colnames(brts_table) <- c("brt", "clade", "event", "endemic") - brts_table[1, ] <- c(1, 0, 0, NA) + brts_table <- matrix(ncol = 5, nrow = 1) + colnames(brts_table) <- c("brt", "clade", "event", "endemic", "col") + brts_table[1, ] <- c(1, 0, 0, NA, NA) expected_IW_format[[1]][[1]] <- list(island_age = 1, not_present = 10, stt_all = stt_all, @@ -254,24 +254,22 @@ test_that("silent with non-empty nonoceanic island with colnames(stt_all) <- c("Time", "nI", "nA", "nC") stt_all[1, ] <- c(1, 1, 2, 0) stt_all[2, ] <- c(0, 0, 2, 0) - brts_table <- matrix(ncol = 4, nrow = 3) - colnames(brts_table) <- c("brt", "clade", "event", "endemic") - brts_table[1, ] <- c(1, 0, 0, NA) - brts_table[2, ] <- c(1, 2, 1, 1) - brts_table[3, ] <- c(1, 1, 1, 1) + brts_table <- matrix(ncol = 5, nrow = 3) + colnames(brts_table) <- c("brt", "clade", "event", "endemic", "col") + brts_table[1, ] <- c(1, 0, 0, NA, NA) + brts_table[2, ] <- c(1, 2, 1, 1, NA) + brts_table[3, ] <- c(1, 1, 1, 1, NA) expected_IW_format[[1]][[1]] <- list(island_age = 1, not_present = 8, stt_all = stt_all, brts_table = brts_table) - expected_IW_format[[1]][[2]] <- list(branching_times = c(1, - 1), + expected_IW_format[[1]][[2]] <- list(branching_times = c(1, 1), stac = 2, - missing_species = 0) - expected_IW_format[[1]][[3]] <- list(branching_times = c(1, - 1), - stac = 2, - missing_species = 0) - expect_true(all.equal(formated_IW_sim, expected_IW_format, tolerance = 1e-7)) + missing_species = 0) + expected_IW_format[[1]][[3]] <- list(branching_times = c(1, 1), + stac = 2, + missing_species = 0) + expect_equal(formated_IW_sim, expected_IW_format, tolerance = 1e-7) }) test_that("silent with non-empty nonoceanic island with @@ -317,52 +315,50 @@ test_that("silent with non-empty nonoceanic island with colnames(stt_all) <- c("Time", "nI", "nA", "nC") stt_all[1, ] <- c(1, 1, 2, 0) stt_all[2, ] <- c(0, 0, 2, 0) - brts_table <- matrix(ncol = 4, nrow = 3) - colnames(brts_table) <- c("brt", "clade", "event", "endemic") - brts_table[1, ] <- c(1, 0, 0, NA) - brts_table[2, ] <- c(1, 2, 1, 1) - brts_table[3, ] <- c(1, 1, 1, 1) + brts_table <- matrix(ncol = 5, nrow = 3) + colnames(brts_table) <- c("brt", "clade", "event", "endemic", "col") + brts_table[1, ] <- c(1, 0, 0, NA, NA) + brts_table[2, ] <- c(1, 2, 1, 1, NA) + brts_table[3, ] <- c(1, 1, 1, 1, NA) expected_IW_format[[1]][[1]] <- list(island_age = 1, not_present = 8, stt_all = stt_all, brts_table = brts_table) - expected_IW_format[[1]][[2]] <- list(branching_times = c(1, - 1), + expected_IW_format[[1]][[2]] <- list(branching_times = c(1, 1), stac = 2, missing_species = 0) - expected_IW_format[[1]][[3]] <- list(branching_times = c(1, - 1), + expected_IW_format[[1]][[3]] <- list(branching_times = c(1, 1), stac = 2, missing_species = 0) - expect_true(all.equal(formated_IW_sim, expected_IW_format, tolerance = 1e-7)) + expect_equal(formated_IW_sim, expected_IW_format, tolerance = 1e-7) }) -test_that("Add_brt_table output is correct when length(island) == 1", { +test_that("add_brt_table output is correct when length(island) == 1", { stt_all <- matrix(ncol = 4, nrow = 2) colnames(stt_all) <- c("Time", "nI", "nA", "nC") stt_all[1, ] <- c(1, 0, 0, 0) stt_all[2, ] <- c(0, 0, 0, 0) island <- list() island[[1]] <- list(island_age = 1, - not_present = 100, - stt_all = stt_all, - init_nonend_spec = 0, - init_end_spec = 0) - formatted_brt <- DAISIE:::Add_brt_table(island) - brt_table <- matrix(ncol = 4, nrow = 1) - colnames(brt_table) <- c("brt", "clade", "event", "endemic") - brt_table[1, ] <- c(1, 0, 0, NA) + not_present = 100, + stt_all = stt_all, + init_nonend_spec = 0, + init_end_spec = 0) + formatted_brt <- DAISIE:::add_brt_table(island) + brt_table <- matrix(ncol = 5, nrow = 1) + colnames(brt_table) <- c("brt", "clade", "event", "endemic", "col") + brt_table[1, ] <- c(1, 0, 0, NA, NA) expected_brt <- list() expected_brt[[1]] <- list(island_age = 1, - not_present = 100, - stt_all = stt_all, - init_nonend_spec = 0, - init_end_spec = 0, - brts_table = brt_table) - expect_true(all.equal(formatted_brt, expected_brt)) + not_present = 100, + stt_all = stt_all, + init_nonend_spec = 0, + init_end_spec = 0, + brts_table = brt_table) + expect_equal(formatted_brt, expected_brt) }) -test_that("Add_brt_table output is correct when length(island) != 1", { +test_that("add_brt_table output is correct when length(island) != 1", { stt_all <- matrix(ncol = 4, nrow = 2) colnames(stt_all) <- c("Time", "nI", "nA", "nC") stt_all[1, ] <- c(1, 0, 0, 0) @@ -387,43 +383,43 @@ test_that("Add_brt_table output is correct when length(island) != 1", { 0.5288428), stac = 4, missing_species = 0) - formatted_brt <- DAISIE:::Add_brt_table(island) - brt_table <- matrix(ncol = 4, nrow = 5) - colnames(brt_table) <- c("brt", "clade", "event", "endemic") - brt_table[1, ] <- c(1, 0, 0, NA) - brt_table[2, ] <- c(0.9244818, 1, 1, 1) - brt_table[3, ] <- c(0.9105857, 1, 2, 1) - brt_table[4, ] <- c(0.5557734, 2, 1, 0) - brt_table[5, ] <- c(0.3146836, 1, 3, 1) + formatted_brt <- DAISIE:::add_brt_table(island) + brt_table <- matrix(ncol = 5, nrow = 5) + colnames(brt_table) <- c("brt", "clade", "event", "endemic", "col") + brt_table[1, ] <- c(1, 0, 0, NA, NA) + brt_table[2, ] <- c(0.9244818, 1, 1, 1, NA) + brt_table[3, ] <- c(0.9105857, 1, 2, 1, NA) + brt_table[4, ] <- c(0.5557734, 2, 1, 0, NA) + brt_table[5, ] <- c(0.3146836, 1, 3, 1, NA) expected_brt <- list() expected_brt[[1]] <- list(island_age = 1, - not_present = 3, - stt_all = stt_all, - init_nonend_spec = 0, - init_end_spec = 0, - brts_table = brt_table) + not_present = 3, + stt_all = stt_all, + init_nonend_spec = 0, + init_end_spec = 0, + brts_table = brt_table) expected_brt[[2]] <- list(branching_times = c(1.0000000, - 0.9244818, - 0.9105857, - 0.3146836), - stac = 2, - missing_species = 0) + 0.9244818, + 0.9105857, + 0.3146836), + stac = 2, + missing_species = 0) expected_brt[[3]] <- list(branching_times = c(1.0000000, - 0.5557734), - stac = 4, - missing_species = 0) + 0.5557734), + stac = 4, + missing_species = 0) expect_equal(formatted_brt, expected_brt) }) -#test_that("Add_brt_table output is correct when length(stac1_5s) != 0") -#test_that("Add_brt_table output is correct when length(stac1_5s) == 0") -#test_that("Add_brt_table output is correct when length(island_no_stac1or5) != 0") +#test_that("add_brt_table output is correct when length(stac1_5s) != 0") +#test_that("add_brt_table output is correct when length(stac1_5s) == 0") +#test_that("add_brt_table output is correct when length(island_no_stac1or5) != 0") test_that("abuse", { expect_error(DAISIE:::DAISIE_format_IW("nonsense")) }) test_that("abuse", { - expect_error(DAISIE:::Add_brt_table("nonsense")) + expect_error(DAISIE:::add_brt_table("nonsense")) }) @@ -477,14 +473,14 @@ test_that("silent with empty island with correct output", { colnames(stt_all) <- c("Time", "nI", "nA", "nC", "nI2", "nA2", "nC2") stt_all[1, ] <- c(1, 0, 0, 0, 0, 0, 0) stt_all[2, ] <- c(0, 0, 0, 0, 0, 0, 0) - brts_table <- matrix(ncol = 4, nrow = 1) - colnames(brts_table) <- c("brt", "clade", "event", "endemic") - brts_table[1, ] <- c(1, 0, 0, NA) + brts_table <- matrix(ncol = 5, nrow = 1) + colnames(brts_table) <- c("brt", "clade", "event", "endemic", "col") + brts_table[1, ] <- c(1, 0, 0, NA, NA) expected_IW_format[[1]][[1]] <- list(island_age = 1, not_present = 20, stt_all = stt_all, brts_table = brts_table) - expect_true(all.equal(formated_IW_sim, expected_IW_format, tolerance = 1e-7)) + expect_equal(formated_IW_sim, expected_IW_format, tolerance = 1e-7) }) test_that("silent when species with two trait states with @@ -543,12 +539,12 @@ test_that("silent when species with two trait states with colnames(stt_all) <- c("Time", "nI", "nA", "nC", "nI2", "nA2", "nC2") stt_all[1, ] <- c(5, 0, 0, 0, 0, 0, 0) stt_all[2, ] <- c(0, 0, 1, 2, 0, 0, 0) - brts_table <- matrix(ncol = 4, nrow = 4) - colnames(brts_table) <- c("brt", "clade", "event", "endemic") - brts_table[1, ] <- c(5.000000000, 0, 0, NA) - brts_table[2, ] <- c(3.102613675, 1, 1, 1) - brts_table[3, ] <- c(1.505629998, 2, 1, 1) - brts_table[4, ] <- c(1.262456559, 2, 2, 1) + brts_table <- matrix(ncol = 5, nrow = 4) + colnames(brts_table) <- c("brt", "clade", "event", "endemic", "col") + brts_table[1, ] <- c(5.000000000, 0, 0, NA, NA) + brts_table[2, ] <- c(3.102613675, 1, 1, 1, NA) + brts_table[3, ] <- c(1.505629998, 2, 1, 1, NA) + brts_table[4, ] <- c(1.262456559, 2, 2, 1, NA) expected_IW_format[[1]][[1]] <- list(island_age = 5, not_present = 13, stt_all = stt_all, diff --git a/tests/testthat/test-DAISIE_loglik_CS.R b/tests/testthat/test-DAISIE_loglik_CS.R index 5660ee2c..dbde56ab 100644 --- a/tests/testthat/test-DAISIE_loglik_CS.R +++ b/tests/testthat/test-DAISIE_loglik_CS.R @@ -13,13 +13,14 @@ test_that("DAISIE_loglik_CS_choice produces correct output for CS_version 1", { brts = brts, stac = stac, missnumspec = missnumspec) + expect_true(is.numeric(loglik)) - expect_equal(loglik, -17.6550433826) + expect_equal(loglik, -17.6535269346579) + }) test_that("DAISIE_loglik_CS_choice produces correct output for relaxed-rate model (CS_version = 2)", { - skip("Produces DLSODES warnings") pars1 <- c(2.000, 2.700, 20.000, 0.009, 1.010) pars2 <- c(1.0e+02, 1.1e+01, 0.0e+00, 0.0e+00, NA, 0.0e+00, 1.0e-04, 1.0e-05, 1.0e-07, 3.0e+03, 9.5e-01, 9.8e-01) @@ -30,51 +31,63 @@ test_that("DAISIE_loglik_CS_choice produces correct output for relaxed-rate CS_version <- list(model = 2, relaxed_par = "cladogenesis", sd = 1) - loglik <- DAISIE_loglik_CS_choice(pars1 = pars1, - pars2 = pars2, - brts = brts, - stac = stac, - missnumspec = missnumspec, - CS_version = CS_version) + + invisible(capture.output(loglik <- DAISIE_loglik_CS_choice(pars1 = pars1, + pars2 = pars2, + brts = brts, + stac = stac, + missnumspec = missnumspec, + CS_version = CS_version))) expect_true(is.numeric(loglik)) - expect_equal(loglik, -9.55117524011) + expect_equal(loglik, -9.550184206825) + }) -test_that("DAISIE_loglik_CS_choice produces correct output for CS_version 0", { +test_that("DAISIE_loglik_CS_choice produces same output for CS_version = 0 (with M = 1) and CS_version = 1 ", { pars1 <- c(2.000, 2.700, 20.000, 0.009, 1.010) - pars2 <- c(1.0e+02, 1.1e+01, 0.0e+00, 0.0e+00, NA, 0.0e+00, 1.0e-04, + pars2 <- c(100, 11, 0, 0, NA, 0.0e+00, 1.0e-04, 1.0e-05, 1.0e-07, 3.0e+03, 9.5e-01, 9.8e-01) brts <- c(4.0000, 3.0282, 1.3227, 0.8223, 0.4286, 0.3462, 0.2450, 0.0808, 0.0527, 0.0327, 0.0221, 0.1180, 0.0756, 0.0525, 0.0322, 0.0118) stac <- 2 missnumspec <- 0 CS_version <- 0 - loglik <- DAISIE_loglik_CS_choice(pars1 = pars1, + datalist <- list(branching_times = brts, stac = stac) + loglik0 <- DAISIE_loglik_CS_choice(pars1 = pars1, pars2 = pars2, + datalist = datalist, brts = brts, stac = stac, missnumspec = missnumspec, CS_version = CS_version) - expect_true(is.numeric(loglik)) - expect_equal(loglik, -17.5608831694) + CS_version <- 1 + loglik1 <- DAISIE_loglik_CS_choice(pars1 = pars1, + pars2 = pars2, + brts = brts, + stac = stac, + missnumspec = missnumspec, + CS_version = CS_version) + + expect_equal(loglik0,loglik1) }) test_that("DAISIE_loglik_all produces correct output for relaxed-rate model", { - skip("Produces DLSODES warnings") utils::data(Galapagos_datalist) - loglik <- DAISIE::DAISIE_loglik_all( - pars1 = c(2.55068735, 2.68345455, 10.00000000, 0.00933207, 1.01007312), - pars2 = c(100, 0, 0, 0, NA), - datalist = Galapagos_datalist, - methode = "lsodes", - CS_version = list(model = 2, - relaxed_par = "cladogenesis", - sd = 1), - abstolint = 1e-16, - reltolint = 1e-10 - ) + invisible(capture.output(suppressWarnings( + loglik <- DAISIE::DAISIE_loglik_all( + pars1 = c(2.55068735, 2.68345455, 10.00000000, 0.00933207, 1.01007312), + pars2 = c(100, 0, 0, 0, NA), + datalist = Galapagos_datalist, + methode = "lsodes", + CS_version = list(model = 2, + relaxed_par = "cladogenesis", + sd = 1), + abstolint = 1e-16, + reltolint = 1e-10 + ) + ))) expect_true(is.numeric(loglik)) - expect_equal(loglik, --77.5108137039949) + expect_equal(loglik, -77.50300644907) }) test_that("DAISIE_loglik produces correct output", { @@ -88,5 +101,49 @@ test_that("DAISIE_loglik produces correct output", { abstolint = 1E-16, reltolint = 1E-10, verbose = FALSE) - expect_equal(output, -0.00347317077256095) + testthat::expect_equal(output, -0.00347317077256095) +}) + +test_that("DAISIE_loglik_all produces same output for CS_version 0 and 1 with and without conditioning", { + utils::data(Galapagos_datalist) + Galapagos_datalist2 <- Galapagos_datalist + for(i in 2:9) { + Galapagos_datalist2[[i]]$branching_times <- c(4, 4 - 2*i*0.1,4 -2*i*0.1-0.1) + Galapagos_datalist2[[i]]$stac <- 2 + } + Galapagos_datalist2 <- DAISIE:::add_brt_table(Galapagos_datalist2) + loglik_CS00 <- DAISIE::DAISIE_loglik_all( + pars1 = c(2.55068735, 2.68345455, 10.00000000, 0.00933207, 1.01007312), + pars2 = c(100, 11, 0, 0, NA), + datalist = Galapagos_datalist2, + methode = "odeint::runge_kutta_fehlberg78", + CS_version = 0, + abstolint = 1e-16, + reltolint = 1e-10) + loglik_CS10 <- DAISIE::DAISIE_loglik_all( + pars1 = c(2.55068735, 2.68345455, 10.00000000, 0.00933207, 1.01007312), + pars2 = c(100, 11, 0, 0, NA), + datalist = Galapagos_datalist2, + methode = "ode45", + CS_version = 1, + abstolint = 1e-16, + reltolint = 1e-10) + testthat::expect_equal(loglik_CS00, loglik_CS10, tol = 5E-6) + loglik_CS01 <- DAISIE::DAISIE_loglik_all( + pars1 = c(2.55068735, 2.68345455, 10.00000000, 0.00933207, 1.01007312), + pars2 = c(100, 11, 1, 0, NA), + datalist = Galapagos_datalist2, + methode = "odeint::runge_kutta_fehlberg78", + CS_version = 0, + abstolint = 1e-16, + reltolint = 1e-10) + loglik_CS11 <- DAISIE::DAISIE_loglik_all( + pars1 = c(2.55068735, 2.68345455, 10.00000000, 0.00933207, 1.01007312), + pars2 = c(100, 11, 1, 0, NA), + datalist = Galapagos_datalist2, + methode = "ode45", + CS_version = 1, + abstolint = 1e-16, + reltolint = 1e-10) + testthat::expect_equal(loglik_CS01, loglik_CS11, tol = 5E-6) }) diff --git a/tests/testthat/test-DAISIE_loglik_integrate.R b/tests/testthat/test-DAISIE_loglik_integrate.R index 9c9c8320..bb60765e 100644 --- a/tests/testthat/test-DAISIE_loglik_integrate.R +++ b/tests/testthat/test-DAISIE_loglik_integrate.R @@ -56,7 +56,7 @@ test_that("DAISIE_loglik_integrate produces correct ouput on radiation", { reltolint = reltolint, verbose = verbose ) - testthat::expect_equal(loglik, -15.1289048939324) + testthat::expect_equal(loglik, -15.12736391328775) }) test_that("DAISIE_loglik_integrand produces correct output", { diff --git a/tests/testthat/test-DAISIE_sim_constant_rate_shift.R b/tests/testthat/test-DAISIE_sim_constant_rate_shift.R index 3b93e73b..effc3673 100644 --- a/tests/testthat/test-DAISIE_sim_constant_rate_shift.R +++ b/tests/testthat/test-DAISIE_sim_constant_rate_shift.R @@ -167,10 +167,12 @@ test_that("Reference output matches DAISIE_sim_constant_rate_shift ", { ) }) -test_that("The SR simulation and inference code works", { - Biwa_datalist <- NULL - rm(Biwa_datalist) - utils::data(Biwa_datalist, package = "DAISIE") +test_that("The SR simulation code works", { + + # Simulate fish diversity over 4 Ma + set.seed(1) + M <- 312 + IslandAge <- 4 pars1 <- c( 0.077, 0.956, @@ -184,16 +186,6 @@ test_that("The SR simulation and inference code works", { 0.442, 0.1951 ) - pars2 <- c(100, 11, 0, 0) - loglik <- DAISIE_SR_loglik_CS(pars1 = pars1, - pars2 = pars2, - datalist = Biwa_datalist) - testthat::expect_equal(loglik, -232.0618, tol = 1E-3) - - # Simulate fish diversity over 4 Ma - set.seed(1) - M <- 312 - IslandAge <- 4 sims <- DAISIE_SR_sim( time = 4, M = M - 17, @@ -206,53 +198,4 @@ test_that("The SR simulation and inference code works", { testthat::expect_equal( unname(sims[[1]][[1]]$stt_all[26, ]), c(0, 56, 11, 0, 66) ) - - Macaronesia_datalist <- NULL - rm(Macaronesia_datalist) - utils::data(Macaronesia_datalist, package = "DAISIE") - pars1 <- c( - 0.1, - 1.1, - 10, - 0.6, - 0.05, - 0.1, - 1.1, - 10, - 0.6, - 0.05, - 7 - ) - pars1mat <- matrix(pars1, nrow = 8, ncol = 11, byrow = T) - expected_loglik <- c( - -Inf, - -266.9125, - -Inf, - -262.7665, - -Inf, - -291.6965, - -Inf, - -261.9222 - ) - # No shift in cladogenesis rate older before - # colonization of diversifying lineages - pars1mat[1, 1] <- 0.01; pars1mat[1, 6] <- 0.3 - pars1mat[2, 1] <- 0.01; pars1mat[2, 6] <- 0.3; pars1mat[2, 11] <- 2.1 - # No shift in extinction rate older than known ages - pars1mat[3, 2] <- 0.1; pars1mat[3, 11] <- 12 - pars1mat[4, 2] <- 0.1; - # No shift in colonization rate older than known ages - pars1mat[5, 4] <- 0.1; pars1mat[5, 9] <- 1.5; pars1mat[5, 11] <- 10 - pars1mat[6, 4] <- 0.1; pars1mat[6, 9] <- 0.8; pars1mat[6, 11] <- 7 - # No shift in anagenetic rate older than any known non-endemic - pars1mat[7, 5] <- 0.01; pars1mat[7, 10] <- 0.1; - pars1mat[8, 5] <- 0.01; pars1mat[8, 10] <- 0.1; pars1mat[8, 11] <- 2.9 - for (i in 1:8) { - loglik <- DAISIE_SR_loglik_CS( - pars1 = pars1mat[i, -12], - pars2 = pars2, - datalist = Macaronesia_datalist[[2]] - ) - testthat::expect_equal(loglik, expected_loglik[i], tol = 1E-3) - } }) diff --git a/tests/testthat/test-integration_DAISIE.R b/tests/testthat/test-integration_DAISIE.R index 54140e13..e46f66c4 100644 --- a/tests/testthat/test-integration_DAISIE.R +++ b/tests/testthat/test-integration_DAISIE.R @@ -19,7 +19,7 @@ test_that("loglik Galapagos works", { ) pars2 <- c(100, 11, 0, 0) loglik <- DAISIE::DAISIE_loglik_all(pars1, pars2, Galapagos_datalist_2types) - testthat::expect_equal(loglik, -61.7094829913735978) + testthat::expect_equal(loglik, -61.70281911731144) }) test_that("loglik macaronesia 2 type works", { @@ -37,7 +37,7 @@ test_that("loglik macaronesia 2 type works", { Macaronesia_datalist[[i]], methode = "lsodes") } - testthat::expect_equal(loglik, -454.9347833283220552) + testthat::expect_equal(loglik, -452.1156003854809) }) test_that("clade specific rate-shift loglik works", { @@ -61,82 +61,68 @@ test_that("clade specific rate-shift loglik works", { }) test_that("IW and CS loglik is same when K = Inf", { - skip_if(Sys.getenv("CI") == "" || !(Sys.getenv("USERNAME") == "rampa"), + skip_if(Sys.getenv("CI") == "" && !(Sys.getenv("USERNAME") == "rampa"), message = "Run only on CI") utils::data(Galapagos_datalist, package = "DAISIE") - pars1 <- c(0.2, 0.1, Inf, 0.001, 0.3) - pars2 <- c(40, 11, 0, 0) - loglik_IW <- DAISIE::DAISIE_loglik_IW( - pars1 = pars1, - pars2 = pars2, - datalist = Galapagos_datalist, - methode = "ode45") - loglik_CS <- DAISIE::DAISIE_loglik_CS( + pars1 <- c(0.35, 0.3, Inf, 0.001, 0.3) + pars2 <- c(80, 11, 1, 0) + Galapagos_datalist_IW <- Galapagos_datalist + for(i in 2:9) { + Galapagos_datalist_IW[[i]]$branching_times <- c(4, 4 - 2*i*0.1,4 -2*i*0.1-0.1) + Galapagos_datalist_IW[[i]]$stac <- 2 + } + Galapagos_datalist_IW <- DAISIE:::add_brt_table(Galapagos_datalist_IW) + invisible(capture.output( + loglik_IW <- DAISIE::DAISIE_loglik_IW( + pars1 = pars1, + pars2 = pars2, + datalist = Galapagos_datalist_IW, + methode = "ode45") + )) + invisible(capture.output( + loglik_IW2 <- DAISIE::DAISIE_loglik_IW( pars1 = pars1, pars2 = pars2, - datalist = Galapagos_datalist, - methode = "ode45", - CS_version = 1) - testthat::expect_lt(abs(loglik_IW - loglik_CS), 5E-6) -}) - -test_that("ontogeny and null-ontogeny loglik is same when ontogeny is - constant", { - skip("Temporary skip") - pars1 <- c(0.2, 0.1, 17, 0.001, 0.3) - pars2 <- c(40, 11, 0, 0) - utils::data(Galapagos_datalist, package = "DAISIE") - loglik_CS <- DAISIE::DAISIE_loglik_all( - pars1 = pars1, - pars2 = pars2, - datalist = Galapagos_datalist, - methode = "ode45") - pars1_td <- c( - max_area = 1, - proportional_peak_t = 0.2, - peak_sharpness = 1, - total_island_age = 15, - lac = pars1[1], - mu_min = pars1[2], - mu_max = pars1[2], - K0 = pars1[3], - gam = pars1[4], - laa = pars1[5] - ) - pars1_td <- DAISIE:::order_pars1(pars1_td) - pars2 <- c(pars2, DAISIE::translate_island_ontogeny("const")) - loglik_time <- DAISIE::DAISIE_loglik_all( - pars1 = pars1_td, - pars2 = pars2, - datalist = Galapagos_datalist, - methode = "ode45" - ) - testthat::expect_equal(loglik_time, loglik_CS) + datalist = Galapagos_datalist_IW, + methode = "odeint::runge_kutta_fehlberg78") + )) + invisible(capture.output( + loglik_CS <- DAISIE::DAISIE_loglik_CS( + pars1 = pars1, + pars2 = pars2, + datalist = Galapagos_datalist_IW, + methode = "ode45", + CS_version = 1) + )) + testthat::expect_equal(loglik_IW, loglik_IW2, tol = 5E-6) + testthat::expect_equal(loglik_IW, loglik_CS, tol = 5E-6) }) testthat::test_that("DAISIE_ML simple case works", { - skip_if(Sys.getenv("CI") == "" || !(Sys.getenv("USERNAME") == "rampa"), + skip_if(Sys.getenv("CI") == "" && !(Sys.getenv("USERNAME") == "rampa"), message = "Run only on CI") expected_mle <- data.frame( - lambda_c = 2.55847849219339, - mu = 2.68768191590176, - K = 6765.0637400135, - gamma = 0.00932987953669849, - lambda_a = 1.00838182578826, - loglik = -76.0001379108545, + lambda_c = 2.583731356303842, + mu = 2.708828027514834, + K = 2992.207701921788, + gamma = 0.00937711049761019, + lambda_a = 0.9993246958280274, + loglik = -75.99266304738612, df = 5L, conv = 0L ) utils::data(Galapagos_datalist) cat("\n") - tested_mle <- DAISIE::DAISIE_ML( - datalist = Galapagos_datalist, - initparsopt = c(2.5, 2.7, 20, 0.009, 1.01), - ddmodel = 11, - idparsopt = 1:5, - parsfix = NULL, - idparsfix = NULL - ) + invisible(capture.output( + tested_mle <- DAISIE::DAISIE_ML( + datalist = Galapagos_datalist, + initparsopt = c(2.5, 2.7, 20, 0.009, 1.01), + ddmodel = 11, + idparsopt = 1:5, + parsfix = NULL, + idparsfix = NULL + ) + )) testthat::expect_equal(expected_mle, tested_mle) }) @@ -147,17 +133,19 @@ test_that("The parameter choice for 2type DAISIE_ML works", { set.seed(1) # MLE and high tolerance for speed-up cat("\n") - fit <- DAISIE::DAISIE_ML( - datalist = Galapagos_datalist_2types, - initparsopt = c(2.183336, 2.517413, 0.009909, 1.080458, 1.316296, 0.001416), - idparsopt = c(1, 2, 4, 5, 7, 11), - parsfix = c(Inf, Inf), - idparsfix = c(3, 8), - idparsnoshift = c(6, 9, 10), - res = 30, - tol = c(1, 1, 1), - maxiter = 30 - ) + invisible(capture.output( + fit <- DAISIE::DAISIE_ML( + datalist = Galapagos_datalist_2types, + initparsopt = c(2.183336, 2.517413, 0.009909, 1.080458, 1.316296, 0.001416), + idparsopt = c(1, 2, 4, 5, 7, 11), + parsfix = c(Inf, Inf), + idparsfix = c(3, 8), + idparsnoshift = c(6, 9, 10), + res = 30, + tol = c(1, 1, 1), + maxiter = 30 + ) + )) testthat::expect_equal(fit$loglik, -74.7557, tol = 1E-3) }) @@ -174,7 +162,7 @@ test_that("conditioning works", { methode = "ode45", CS_version = 1 ) - testthat::expect_equal(loglik_CS_1type_cond0, -96.49629968062564) + testthat::expect_equal(loglik_CS_1type_cond0, -96.49069330275196) ## 2 type utils::data(Galapagos_datalist_2types, package = "DAISIE") @@ -197,7 +185,7 @@ test_that("conditioning works", { pars2_2type_cond0, Galapagos_datalist_2types ) - testthat::expect_equal(loglik_CS_2type_cond0, -61.709482984890265) + testthat::expect_equal(loglik_CS_2type_cond0, -61.70281911731144) # Cond 1 ## 1 type @@ -211,7 +199,7 @@ test_that("conditioning works", { methode = 'ode45', CS_version = 1 ) - testthat::expect_equal(loglik_CS_1type_cond1, -96.463184608046333) + testthat::expect_equal(loglik_CS_1type_cond1, -96.45757823017264) ## 2 type utils::data(Galapagos_datalist_2types, package = "DAISIE") @@ -234,7 +222,7 @@ test_that("conditioning works", { pars2_2type_cond1, Galapagos_datalist_2types ) - testthat::expect_equal(loglik_CS_2type_cond1,-61.4442595468189054) + testthat::expect_equal(loglik_CS_2type_cond1, -61.4375956792401) # Cond 5 ## 1 type @@ -248,7 +236,7 @@ test_that("conditioning works", { methode = 'ode45', CS_version = 1 ) - testthat::expect_equal(loglik_CS_1type_cond5,-95.1456087499799423) + testthat::expect_equal(loglik_CS_1type_cond5, -95.14000237210625) ## 2 type utils::data(Galapagos_datalist_2types, package = "DAISIE") @@ -271,5 +259,5 @@ test_that("conditioning works", { pars2_2type_cond5, Galapagos_datalist_2types ) - testthat::expect_equal(loglik_CS_2type_cond5, -61.3801835140081025) + testthat::expect_equal(loglik_CS_2type_cond5, -61.3735196464293) }) diff --git a/vignettes/demo_optimize.Rmd b/vignettes/demo_optimize.Rmd index daf15d97..dc228050 100644 --- a/vignettes/demo_optimize.Rmd +++ b/vignettes/demo_optimize.Rmd @@ -22,7 +22,7 @@ Citation: Valente LM, Phillimore AB, Etienne RS (2015) Equilibrium and non- equilibrium dynamics simultaneously operate in the Galápagos islands. Ecology Letters, In press. -## 3.2 Data table +## Loading data table To load the package: @@ -75,7 +75,7 @@ The same data can also be visualized: DAISIE::DAISIE_plot_island(Galapagos_datatable, island_age = 4) ``` -## 3.3 Formatting table to run in DAISIE +## Formatting table to run in DAISIE Before running analyses, the datatable needs to be converted to a DAISIE datalist format using the function DAISIE_dataprep. @@ -117,7 +117,7 @@ Galapagos_datalist_2types <- DAISIE_dataprep( The objects `Galapagos_datalist` and `Galapagos_datalist_2types` can now be run directly in `DAISIE` functions. -## 3.4 Optimizing parameters using maximum likelihood +## Optimizing parameters using maximum likelihood The function that conducts maximum likelihood optimization of DAISIE model parameters is diff --git a/vignettes/demo_optimize.html b/vignettes/demo_optimize.html index aee6d9e2..59bd3408 100644 --- a/vignettes/demo_optimize.html +++ b/vignettes/demo_optimize.html @@ -16,23 +16,40 @@ Demo optimizing parameters + + - + + + + @@ -310,19 +148,19 @@

19 May 2015

DAISIE – Dynamic Assembly of Island biotas through Speciation, Immigration and Extinction

Citation: Valente LM, Phillimore AB, Etienne RS (2015) Equilibrium and non- equilibrium dynamics simultaneously operate in the Galápagos islands. Ecology Letters, In press.

-
-

3.2 Data table

+
+

Loading data table

To load the package:

-
library(DAISIE)
+
library(DAISIE)

The raw dataset is inputted as a table. The Galápagos dataset table can be visualized with:

-
data(Galapagos_datatable)
-knitr::kable(Galapagos_datatable)
+
data(Galapagos_datatable)
+knitr::kable(Galapagos_datatable)
----++++ @@ -397,41 +235,41 @@

3.2 Data table

  • Branching_times – This should be the 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 branching times of the radiation including the stem age of the radiation. Note – if there are species within the radiation that are not found on the island (e.g. back-colonisation) the branching times of these species should be excluded, as the mainland species pool is treated as static.

    • The same data can also be visualized:

    -
    DAISIE::DAISIE_plot_island(Galapagos_datatable, island_age = 4)
-#> [1] "Colonisation time of 7.456 for Coccyzus is older than island age"
-#> [1] "Colonisation time of 10.285 for Pyrocephalus is older than island age"
    +
    DAISIE::DAISIE_plot_island(Galapagos_datatable, island_age = 4)
+#> [1] "Colonisation time of 7.456 for Coccyzus is older than island age"
+#> [1] "Colonisation time of 10.285 for Pyrocephalus is older than island age"

    -

    3.3 Formatting table to run in DAISIE

    +

    Formatting table to run in DAISIE

    Before running analyses, the datatable needs to be converted to a DAISIE datalist format using the function DAISIE_dataprep.

    We will prepare two different datalists based on the Galápagos datatable. In the 1st datalist we will treat all taxa as equivalent. We will specify an island age of four million years (island_age=4) and a mainland pool size of 1000 (M=1000).

    -
    data(Galapagos_datatable) 
-      
-Galapagos_datalist <- DAISIE_dataprep( 
-  datatable = Galapagos_datatable, 
-  island_age = 4, 
-  M = 1000
-)
-#> [1] "Colonisation time of 7.456 for Coccyzus is older than island age"
-#> [1] "Colonisation time of 10.285 for Pyrocephalus is older than island age"
    +
    data(Galapagos_datatable) 
+      
+Galapagos_datalist <- DAISIE_dataprep( 
+  datatable = Galapagos_datatable, 
+  island_age = 4, 
+  M = 1000
+)
+#> [1] "Colonisation time of 7.456 for Coccyzus is older than island age"
+#> [1] "Colonisation time of 10.285 for Pyrocephalus is older than island age"

    In the 2nd datalist we will allow for the Darwin’s finches to form a separate group for which rates can be decoupled from those governing the macroevolutionary process in all other clades (number_clade_types=2 and list_type2_clades = “Finches”). We will set the proportion of Darwin’s finch type species in the mainland pool to be 0.163. (prop_type2_pool=0.163). If prop_type2_pool is not specified then by default it is given the value of the proportion of the Galapagos lineages that Darwin’s finches represent (1/8=0.125 in this case).

    -
    data(Galapagos_datatable) 
-
-Galapagos_datalist_2types <- DAISIE_dataprep( 
-  datatable = Galapagos_datatable, 
-  island_age = 4, 
-  M = 1000, 
-  number_clade_types = 2, 
-  list_type2_clades = "Finches", 
-  prop_type2_pool = 0.163
-)
-#> [1] "Colonisation time of 7.456 for Coccyzus is older than island age"
-#> [1] "Colonisation time of 10.285 for Pyrocephalus is older than island age"
    +
    data(Galapagos_datatable) 
+
+Galapagos_datalist_2types <- DAISIE_dataprep( 
+  datatable = Galapagos_datatable, 
+  island_age = 4, 
+  M = 1000, 
+  number_clade_types = 2, 
+  list_type2_clades = "Finches", 
+  prop_type2_pool = 0.163
+)
+#> [1] "Colonisation time of 7.456 for Coccyzus is older than island age"
+#> [1] "Colonisation time of 10.285 for Pyrocephalus is older than island age"

    The objects Galapagos_datalist and Galapagos_datalist_2types can now be run directly in DAISIE functions.

    -

    3.4 Optimizing parameters using maximum likelihood

    +

    Optimizing parameters using maximum likelihood

    The function that conducts maximum likelihood optimization of DAISIE model parameters is called DAISIE_ML.

    Different models can be specified using ddmodel option in DAISIE_ML:

      diff --git a/vignettes/demo_sim.Rmd b/vignettes/demo_sim.Rmd index 7f34b21c..5bf9026b 100644 --- a/vignettes/demo_sim.Rmd +++ b/vignettes/demo_sim.Rmd @@ -22,7 +22,7 @@ Citation: Valente LM, Phillimore AB, Etienne RS (2015) Equilibrium and non- equilibrium dynamics simultaneously operate in the Galápagos islands. Ecology Letters, 18(8), 844-852. https://doi.org/10.1111/ele.12461. -## 3.2 Data table +## Loading To load the package: @@ -30,7 +30,7 @@ To load the package: library(DAISIE) ``` -## 3.5 Simulating islands +## Simulating islands The function `DAISIE_sim` allows simulation of DAISIE models and plots the results. The user specifies the parameters to be simulated, the number of replicates, the length of the diff --git a/vignettes/demo_sim.html b/vignettes/demo_sim.html index 265dba34..66b826ed 100644 --- a/vignettes/demo_sim.html +++ b/vignettes/demo_sim.html @@ -16,23 +16,40 @@ Demo simulating islands + + - + + + + @@ -310,87 +148,87 @@

      19 May 2015

      DAISIE – Dynamic Assembly of Island biotas through Speciation, Immigration and Extinction

      Citation: Valente LM, Phillimore AB, Etienne RS (2015) Equilibrium and non- equilibrium dynamics simultaneously operate in the Galápagos islands. Ecology Letters, 18(8), 844-852. https://doi.org/10.1111/ele.12461.

      -
      -

      3.2 Data table

      +
      +

      Loading

      To load the package:

      -
      library(DAISIE)
      +
      library(DAISIE)
      -

      3.5 Simulating islands

      +

      Simulating islands

      The function DAISIE_sim allows simulation of DAISIE models and plots the results. The user specifies the parameters to be simulated, the number of replicates, the length of the simulation (typically the island age), and the number of species in the mainland pool.

      When the plot_sims option is set to the default (TRUE) the function will produce a species- through-time plot showing the accumulation of total, endemic and non-endemic species through time, as well as confidence intervals for the total number of species.

      -
      n_mainland_species <- 1000
-island_age <- 4
      +
      n_mainland_species <- 1000
+island_age <- 4

      To shorten the run-time of this vignette, reduce the number of n_replicates. For increased accuracy, increase this number.

      -
      n_replicates <- 10
      +
      n_replicates <- 10

      Example 5.1 – Simulating 10 islands with no diversity-dependence, all species sharing the same parameters, and plotting the results

      -
      set.seed(42)
      -
      clado_rate <- 2.550687345 # cladogenesis rate
-ext_rate <- 2.683454548 # extinction rate
-clade_carr_cap <- Inf # clade-level carrying capacity
-imm_rate <- 0.00933207 # immigration rate
-ana_rate <- 1.010073119 # anagenesis rate
-
-island_replicates <- DAISIE_sim_constant_rate( 
-  time = island_age,
-  M = n_mainland_species, 
-  pars = c(clado_rate, ext_rate, clade_carr_cap, imm_rate, ana_rate),
-  replicates = n_replicates,
-  plot_sims = FALSE,
-  verbose = FALSE
-)
-DAISIE_plot_sims(island_replicates = island_replicates)
      +
      set.seed(42)
      +
      clado_rate <- 2.550687345 # cladogenesis rate
+ext_rate <- 2.683454548 # extinction rate
+clade_carr_cap <- Inf # clade-level carrying capacity
+imm_rate <- 0.00933207 # immigration rate
+ana_rate <- 1.010073119 # anagenesis rate
+
+island_replicates <- DAISIE_sim_constant_rate( 
+  time = island_age,
+  M = n_mainland_species, 
+  pars = c(clado_rate, ext_rate, clade_carr_cap, imm_rate, ana_rate),
+  replicates = n_replicates,
+  plot_sims = FALSE,
+  verbose = FALSE
+)
+DAISIE_plot_sims(island_replicates = island_replicates)

      The object island_replicates contains the results of the simulation in DAISIE format. 10 islands are stored in the object, and each island replicate can be viewed separately. For example, type island_replicates[[1]] to view the first replicate.

      The element of the list relating to each island contains a table with the number of species through time, as well as branching time information for each independent colonisation event extant at the end of the simulation.

      Example 5.2 – Simulating 10 islands with diversity-dependence (K’=10), all species sharing the same parameters, and plotting the results

      -
      clado_rate <- 2.550687345 # cladogenesis rate
-ext_rate <- 2.683454548 # extinction rate
-clade_carr_cap <- 10.0  # clade-level carrying capacity
-imm_rate <- 0.00933207 # immigration rate
-ana_rate <- 1.010073119 # anagenesis rate
-
-island_replicates_K <- DAISIE_sim_constant_rate( 
-  time = island_age, 
-  M = n_mainland_species, 
-  pars = c(clado_rate, ext_rate, clade_carr_cap, imm_rate, ana_rate),
-  replicates = n_replicates,
-  plot_sims = FALSE,
-  verbose = FALSE
-) 
-DAISIE_plot_sims(island_replicates_K)
      +
      clado_rate <- 2.550687345 # cladogenesis rate
+ext_rate <- 2.683454548 # extinction rate
+clade_carr_cap <- 10.0  # clade-level carrying capacity
+imm_rate <- 0.00933207 # immigration rate
+ana_rate <- 1.010073119 # anagenesis rate
+
+island_replicates_K <- DAISIE_sim_constant_rate( 
+  time = island_age, 
+  M = n_mainland_species, 
+  pars = c(clado_rate, ext_rate, clade_carr_cap, imm_rate, ana_rate),
+  replicates = n_replicates,
+  plot_sims = FALSE,
+  verbose = FALSE
+) 
+DAISIE_plot_sims(island_replicates_K)

      Example 5.3 – Simulating 10 islands allowing Darwin’s finches to have a higher rate of cladogenesis:

      -
      clado_rate_1 <- 0.38 # cladogenesis rate
-ext_rate_1 <- 0.55 # extinction rate
-clade_carr_cap_1 <- Inf  # clade-level carrying capacity
-imm_rate_1 <- 0.04 # immigration rate
-ana_rate_1 <- 1.10 # anagenesis rate
-
-clado_rate_2 <- 0.38 # cladogenesis rate
-ext_rate_2 <- ext_rate_1 # extinction rate
-clade_carr_cap_2 <- clade_carr_cap_1  # clade-level carrying capacity
-imm_rate_2 <- imm_rate_1 # immigration rate
-ana_rate_2 <- ana_rate_1 # anagenesis rate
-
-island_replicates_2types <- DAISIE_sim_constant_rate( 
-  time = island_age,
-  M = n_mainland_species, 
-  pars = c(
-    clado_rate_1, ext_rate_1, clade_carr_cap_1, imm_rate_1, ana_rate_1,
-    clado_rate_2, ext_rate_2, clade_carr_cap_2, imm_rate_2, ana_rate_2
-  ),
-  replicates = n_replicates, 
-  prop_type2_pool = 0.163,
-  plot_sims = FALSE,
-  verbose = FALSE
-)
-DAISIE_plot_sims(island_replicates_2types)
      +
      clado_rate_1 <- 0.38 # cladogenesis rate
+ext_rate_1 <- 0.55 # extinction rate
+clade_carr_cap_1 <- Inf  # clade-level carrying capacity
+imm_rate_1 <- 0.04 # immigration rate
+ana_rate_1 <- 1.10 # anagenesis rate
+
+clado_rate_2 <- 0.38 # cladogenesis rate
+ext_rate_2 <- ext_rate_1 # extinction rate
+clade_carr_cap_2 <- clade_carr_cap_1  # clade-level carrying capacity
+imm_rate_2 <- imm_rate_1 # immigration rate
+ana_rate_2 <- ana_rate_1 # anagenesis rate
+
+island_replicates_2types <- DAISIE_sim_constant_rate( 
+  time = island_age,
+  M = n_mainland_species, 
+  pars = c(
+    clado_rate_1, ext_rate_1, clade_carr_cap_1, imm_rate_1, ana_rate_1,
+    clado_rate_2, ext_rate_2, clade_carr_cap_2, imm_rate_2, ana_rate_2
+  ),
+  replicates = n_replicates, 
+  prop_type2_pool = 0.163,
+  plot_sims = FALSE,
+  verbose = FALSE
+)
+DAISIE_plot_sims(island_replicates_2types)

      This produces a figure similar to Fig. 2, with three plots: one for the total number of species, one for species of type 1 and one for species of type 2. Accessing each island replicate individually (e.g. island_replicates_2types[[15]]) shows information on branching times and species-through-time tables for total, type 1 species and type 2 species.