diff --git a/.Rbuildignore b/.Rbuildignore index d8c1274a..5d79eddc 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -15,3 +15,4 @@ ^\.Rproj\.user$ ^\.covrignore$ ^LICENSE$ +^\.zenodo\.json$ diff --git a/.zenodo.json b/.zenodo.json new file mode 100644 index 00000000..67221f93 --- /dev/null +++ b/.zenodo.json @@ -0,0 +1,104 @@ +{ + "creators":[ + { + "orcid":"0000-0003-2142-7612", + "affiliation":"Groningen Institute for Evolutionary Life Sciences, University of Groningen", + "name":"Etienne, Rampal S." + }, + { + "orcid":"0000-0003-4247-8785", + "affiliation":"Naturalis Biodiversity Center", + "name":"Valente, Luis" + }, + { + "orcid":"0000-0002-6553-1553", + "affiliation":"Institute of Evolutionary Biology, University of Edinburgh", + "name":"Phillimore, Albert B." + }, + { + "orcid":"0000-0003-2325-4727", + "affiliation":"Station d’Ecologie Théorique et Expérimentale, CNRS", + "name":"Haegeman, Bart" + }, + { + "orcid":"0000-0001-5218-3046", + "affiliation":"Groningen Institute for Evolutionary Life Sciences, University of Groningen", + "name":"Lambert, Joshua W." + }, + { + "orcid":"0000-0003-2561-4677", + "affiliation":"Groningen Institute for Evolutionary Life Sciences, University of Groningen", + "name":"Santos Neves, Pedro" + }, + { + "orcid":"0000-0001-9594-946X", + "affiliation":"Groningen Institute for Evolutionary Life Sciences, University of Groningen", + "name":"Xie, Shu" + }, + { + "orcid":"0000-0003-1107-7049", + "affiliation":"Groningen Institute for Evolutionary Life Sciences, University of Groningen", + "name":"Bilderbeek, Rich\u00e8l J.C." + }, + { + "orcid":"0000-0002-6784-1037", + "affiliation":"Groningen Institute for Evolutionary Life Sciences, University of Groningen", + "name":"Hildenbrandt, Hanno" + } + ], + "contributors":[ + { + "orcid":"0000-0001-5711-9457", + "affiliation":"Justus Liebig University", + "name":"Hauffe, Torsten", + "type":"Other" + }, + { + "orcid":"0000-0002-2952-3345", + "affiliation":"Groningen Institute for Evolutionary Life Sciences, University of Groningen", + "name":"Laudanno, Giovanni", + "type":"Other" + }, + { + "orcid":"0000-0002-9720-4581", + "affiliation":"National University of Singapore", + "name":"Kristensen, Nadiah", + "type":"Other" + }, + { + "orcid":"0000-0002-1447-7630", + "affiliation":"Groningen Institute for Evolutionary Life Sciences, University of Groningen", + "name":"Scherrer, Raphael", + "type":"Other" + } + ], + "license":"GPL-3.0-only", + "title":"DAISIE: Dynamical Assembly of Islands by Speciation, Immigration and Extinction", + "description":"Simulates and computes the (maximum) likelihood of a dynamical model of island biota assembly through speciation, immigration and extinction.", + "access_right":"open", + "keywords":[ + "island biogeography", + "birth-death process", + "maximum-likelihood estimation", + "R package" + ], + "related_identifiers":[ + { + "scheme":"url", + "identifier":"https://github.com/rsetienne/DAISIE", + "relation":"isSupplementTo", + "resource_type":"software" + }, + { + "scheme":"url", + "identifier":"https://cran.r-project.org/package=DAISIE", + "relation":"isSupplementTo", + "resource_type":"software" + }, + { + "scheme":"doi", + "identifier":"10.5281/zenodo.4054058", + "relation":"isVersionOf" + } + ] +} \ No newline at end of file diff --git a/DESCRIPTION b/DESCRIPTION index 8e6cf02e..00b1d0fa 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,8 +1,8 @@ Package: DAISIE Type: Package Title: Dynamical Assembly of Islands by Speciation, Immigration and Extinction -Version: 4.1.1 -Date: 2022-04-14 +Version: 4.2.0 +Date: 2022-05-24 Depends: R (>= 3.5.0) biocViews: SystemRequirements: C++14 @@ -115,5 +115,4 @@ Encoding: UTF-8 VignetteBuilder: knitr URL: https://github.com/rsetienne/DAISIE BugReports: https://github.com/rsetienne/DAISIE/issues -RoxygenNote: 7.1.2 -LazyData: true +RoxygenNote: 7.2.0 diff --git a/NAMESPACE b/NAMESPACE index 40da16c6..49285890 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -12,6 +12,7 @@ export(DAISIE_SR_ML) export(DAISIE_SR_ML_CS) export(DAISIE_SR_loglik_CS) export(DAISIE_SR_loglik_all) +export(DAISIE_abm_factor) export(DAISIE_convertprobdist) export(DAISIE_dataprep) export(DAISIE_loglik_CS) diff --git a/NEWS.md b/NEWS.md index 8fbf6be0..27ac4fbc 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,24 @@ +# DAISIE 4.2.0 + +* Important bugfixes on estimation when data contains a lineage or a clade for +which (see `stac_key` vignette for details on each case): + * The divergence time is unknown but known to have occurred *after* a + specific point in time in the island's existence. + * The divergence time is unknown but known to have occurred *before* a + specific point in time in the island's existence (usually known when dated + population level phylogenetic data is available). +* Improve and finalise the continental island estimation scenario, when there +are initially species present on the island on inception. +* Improve tests related to the two points above. +* Remove DAISIE time dependent estimation code for now to ease work on the remaining +code. Development will proceed separately. +* Rate computations in simulation are first calculated per lineage, to then be +calculated per island. This will be needed for time dependent estimation later on. +* Improve zenodo metadata. +* Add `methode = "odeint::adams_bashforth_moulton_X"` to the list of available +numeric integrators for parameter estimation. + + # DAISIE 4.1.1 * Correctly use `is.data.frame()` rather than `class(foo) == "data.frame"` to satisfy CRAN note. diff --git a/R/DAISIE_ML1.R b/R/DAISIE_ML1.R index 6b0dad08..9eacd329 100644 --- a/R/DAISIE_ML1.R +++ b/R/DAISIE_ML1.R @@ -11,18 +11,28 @@ DAISIE_loglik_all_choosepar <- function(trparsopt, CS_version = 1, abstolint = 1E-16, reltolint = 1E-10) { - if (sum(idparsnoshift %in% (6:10)) != 5) { + all_no_shift <- 6:10 + if (max(idparsopt,-Inf) <= 6 && + max(idparsfix,-Inf) <= 6 && + (6 %in% idparsopt || 6 %in% idparsfix)) { + idparsnoshift <- 7:11 + all_no_shift <- 7:11 + } + if (sum(idparsnoshift %in% (all_no_shift)) != 5) { trpars1 <- rep(0, 11) } else { - trpars1 <- rep(0, 5) - trparsfix <- trparsfix[-which(idparsfix == 11)] - idparsfix <- idparsfix[-which(idparsfix == 11)] + trpars1 <- rep(0, 6) + prop_type2_present <- which(idparsfix == 11) + if(!is.null(prop_type2_present)) { + trparsfix <- trparsfix[-prop_type2_present] + idparsfix <- idparsfix[-prop_type2_present] + } } trpars1[idparsopt] <- trparsopt if (length(idparsfix) != 0) { trpars1[idparsfix] <- trparsfix } - if (sum(idparsnoshift %in% (6:10)) != 5) { + if (sum(idparsnoshift %in% all_no_shift) != 5) { trpars1[idparsnoshift] <- trpars1[idparsnoshift - 5] } if (max(trpars1) > 1 | min(trpars1) < 0) { @@ -31,14 +41,14 @@ DAISIE_loglik_all_choosepar <- function(trparsopt, pars1 <- trpars1 / (1 - trpars1) if (pars2[6] > 0) { pars1 <- DAISIE_eq(datalist, pars1, pars2[-5]) - if (sum(idparsnoshift %in% (6:10)) != 5) { + if (sum(idparsnoshift %in% all_no_shift) != 5) { pars1[idparsnoshift] <- pars1[idparsnoshift - 5] } } if (min(pars1) < 0) { loglik <- -Inf } else { - loglik <- DAISIE::DAISIE_loglik_all( + loglik <- DAISIE_loglik_all( pars1 = pars1, pars2 = pars2, datalist = datalist, @@ -180,7 +190,35 @@ DAISIE_ML1 <- function( "lambda_a2", "prop_type2" ) + all_no_shift <- 6:10 + max_idpars <- 11 + if (max(idparsopt, -Inf) <= 6 && + max(idparsfix, -Inf) <= 6 && + (6 %in% idparsopt || 6 %in% idparsfix)) { + max_idpars <- 12 + idparsnoshift <- 7:11 + all_no_shift <- 7:11 + namepars <- c( + "lambda_c", + "mu", + "K", + "gamma", + "lambda_a", + "prob_init_pres", + "lambda_c2", + "mu2", + "K2", + "gamma2", + "lambda_a2", + "prop_type2" + ) + nc <- NA + names(nc) <- "prob_init_pres" + out2err <- add_column_to_dataframe(df = out2err, + position = 'lambda_a', + column_to_insert = nc) + } if (length(namepars[idparsopt]) == 0) { optstr <- "nothing" } else { @@ -194,16 +232,13 @@ DAISIE_ML1 <- function( fixstr <- namepars[idparsfix] } cat("You are fixing", fixstr, "\n") - if (sum(idparsnoshift %in% (6:10)) != 5) { + + if (sum(idparsnoshift %in% (all_no_shift)) != 5) { noshiftstring <- namepars[idparsnoshift] cat("You are not shifting", noshiftstring, "\n") } idpars <- sort(c(idparsopt, idparsfix, idparsnoshift, idparseq)) - if (!any(idpars == 11)) { - idpars <- c(idpars, 11) - idparsfix <- c(idparsfix, 11) - parsfix <- c(parsfix, 0) - } + missnumspec <- unlist(lapply(datalist, function(list) {list$missing_species})) # nolint if (sum(missnumspec) > (res - 1)) { cat( @@ -211,11 +246,13 @@ DAISIE_ML1 <- function( resolution of the ODE.\n") return(out2err) } - if (length(idpars) != 11) { - cat("You have too many parameters to be optimized or fixed.\n") + + if ((length(idpars) != max(idpars))) { + cat("The parameters to be optimized and/or fixed are incoherent.\n") return(out2err) } - if ((prod(idpars == (1:11)) != 1) || # nolint + + if ((!all(idpars == 1:max(idpars))) || # nolint (length(initparsopt) != length(idparsopt)) || (length(parsfix) != length(idparsfix))) { cat("The parameters to be optimized and/or fixed are incoherent.\n") @@ -313,10 +350,10 @@ DAISIE_ML1 <- function( MLpars <- MLtrpars / (1 - MLtrpars) ML <- as.numeric(unlist(out$fvalues)) - if (sum(idparsnoshift %in% (6:10)) != 5) { - MLpars1 <- rep(0, 10) + if (sum(idparsnoshift %in% (all_no_shift)) != 5) { + MLpars1 <- rep(0, 11) } else { - MLpars1 <- rep(0, 5) + MLpars1 <- rep(0, 6) } MLpars1[idparsopt] <- MLpars if (length(idparsfix) != 0) { @@ -331,7 +368,7 @@ DAISIE_ML1 <- function( MLpars1[3] <- Inf } - if (sum(idparsnoshift %in% (6:10)) != 5) { + if (sum(idparsnoshift %in% (all_no_shift)) != 5) { if (length(idparsnoshift) != 0) { MLpars1[idparsnoshift] <- MLpars1[idparsnoshift - 5] } @@ -355,9 +392,7 @@ DAISIE_ML1 <- function( conv = unlist(out$conv) ) s1 <- sprintf( - "Maximum likelihood parameter estimates: lambda_c: %f, mu: %f, K: %f, - gamma: %f, lambda_a: %f, lambda_c2: %f, mu2: %f, K2: %f, gamma2: %f, - lambda_a2: %f, prop_type2: %f", + "Maximum likelihood parameter estimates:\n lambda_c: %f\n mu: %f\n K: %f\n gamma: %f\n lambda_a: %f\n lambda_c2: %f\n mu2: %f\n K2: %f\n gamma2: %f\n lambda_a2: %f\n prop_type2: %f", MLpars1[1], MLpars1[2], MLpars1[3], @@ -370,6 +405,27 @@ DAISIE_ML1 <- function( MLpars1[10], MLpars1[11] ) + } else if (all(all_no_shift == 7:11)) { + out2 <- data.frame( + lambda_c = MLpars1[1], + mu = MLpars1[2], + K = MLpars1[3], + gamma = MLpars1[4], + lambda_a = MLpars1[5], + prob_init_pres = MLpars1[6], + loglik = ML, + df = length(initparsopt), + conv = unlist(out$conv) + ) + s1 <- sprintf( + "Maximum likelihood parameter estimates:\n lambda_c: %f\n mu: %f\n K: %f\n gamma: %f\n lambda_a: %f\n prob_init_pres: %f", + MLpars1[1], + MLpars1[2], + MLpars1[3], + MLpars1[4], + MLpars1[5], + MLpars1[6] + ) } else { out2 <- data.frame( lambda_c = MLpars1[1], @@ -382,8 +438,7 @@ DAISIE_ML1 <- function( conv = unlist(out$conv) ) s1 <- sprintf( - "Maximum likelihood parameter estimates: lambda_c: %f, mu: %f, K: %f, - gamma: %f, lambda_a: %f", + "Maximum likelihood parameter estimates:\n lambda_c: %f\n mu: %f\n K: %f\n gamma: %f\n lambda_a: %f\n", MLpars1[1], MLpars1[2], MLpars1[3], diff --git a/R/DAISIE_ML3.R b/R/DAISIE_ML3.R deleted file mode 100644 index 0017b46c..00000000 --- a/R/DAISIE_ML3.R +++ /dev/null @@ -1,191 +0,0 @@ -# Don't document this function. For internal use only. -DAISIE_loglik_all_choosepar3 = function( - trparsopt, - trparsfix, - idparsopt, - idparsfix, - pars2, - datalist, - methode, - CS_version = 1, - abstolint = 1E-16, - reltolint = 1E-10 - ) { - trpars1 <- rep(0, 10) - trpars1[idparsopt] <- trparsopt - if (length(idparsfix) != 0) { - trpars1[idparsfix] <- trparsfix - } - if (max(trpars1) > 1 | min(trpars1) < 0) { - loglik <- -Inf - } else { - pars1 <- trpars1 / (1 - trpars1) - if (min(pars1) < 0) { - loglik <- -Inf - } else { - loglik <- DAISIE_loglik_all(pars1 = pars1, pars2 = pars2, datalist = datalist, methode = methode, CS_version = CS_version, abstolint = abstolint, reltolint = reltolint) - } - if (is.nan(loglik) || is.na(loglik)) { - cat("There are parameter values used which cause numerical problems.\n") - loglik <- -Inf - } - } - return(loglik) -} - -#' Computes MLE for single type species under a clade specific scenario with -#' ontogeny -#' -#' @inheritParams default_params_doc -#' -#' @return The output is a dataframe containing estimated parameters and -#' maximum loglikelihood. \item{lambda_c}{ gives the maximum likelihood -#' estimate of lambda^c, the rate of cladogenesis} \item{mu}{ gives the maximum -#' likelihood estimate of mu, the extinction rate} \item{K}{ gives the maximum -#' likelihood estimate of K, the carrying-capacity} \item{gamma}{ gives the -#' maximum likelihood estimate of gamma, the immigration rate } -#' \item{lambda_a}{ gives the maximum likelihood estimate of lambda^a, the rate -#' of anagenesis} \item{loglik}{ gives the maximum loglikelihood} -#' \item{df}{ gives the number of estimated parameters, i.e. degrees of feedom} -#' \item{conv}{ gives a message on convergence of optimization; -#' conv = 0 means convergence} -#' @keywords internal -DAISIE_ML3 <- function( - datalist, - initparsopt, - idparsopt, - parsfix, - idparsfix, - res = 100, - ddmodel = 0, - cond = 0, - island_ontogeny, - tol = c(1E-4, 1E-5, 1E-7), - maxiter = 1000 * round((1.25) ^ length(idparsopt)), - methode = "lsodes", - optimmethod = "subplex", - CS_version = 1, - verbose = 0, - tolint = c(1E-16, 1E-10), - 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 -# - next largest brts = stem age / time of divergence from the mainland -# The interpretation of this depends on stac (see below) -# For stac = 0, this needs to be specified only once. -# For stac = 1, this is the time since divergence from the immigrant's sister on the mainland. -# The immigrant must have immigrated at some point since then. -# For stac = 2 and stac = 3, this is the time since divergence from the mainland. -# The immigrant that established the clade on the island must have immigrated precisely at this point. -# For stac = 3, it must have reimmigrated, but only after the first immigrant had undergone speciation. -# - min(brts) = most recent branching time (only for stac = 2, or stac = 3) -# datalist[[,]][2] = list of status of the clades formed by the immigrant -# . stac == 0 : immigrant is not present and has not formed an extant clade -# Instead of a list of zeros, here a number must be given with the number of clades having stac = 0 -# . stac == 1 : immigrant is present but has not formed an extant clade -# . stac == 2 : immigrant is not present but has formed an extant clade -# . stac == 3 : immigrant is present and has formed an extant clade -# . stac == 4 : immigrant is present but has not formed an extant clade, and it is known when it immigrated. -# . stac == 5 : immigrant is not present and has not formed an extent clade, but only an endemic species -# datalist[[,]][3] = list with number of missing species in clades for stac = 2 and stac = 3; -# for stac = 0 and stac = 1, this number equals 0. -# initparsopt, parsfix = optimized and fixed model parameters -# - pars1[1:4] = area_pars -# - pars1[5] = lac = (initial) cladogenesis rate -# - pars1[6:7] = extinction rate parameters -# - pars1[8] = K = maximum number of species possible in the clade -# - pars1[9] = gam = (initial) immigration rate -# - pars1[10] = laa = (initial) anagenesis rate -# idparsopt, idparsfix = ids of optimized and fixed model parameters -# - res = pars2[1] = lx = length of ODE variable x -# - ddmodel = pars2[2] = diversity-dependent model,mode of diversity-dependence -# . ddmodel == 0 : no diversity-dependence -# . ddmodel == 1 : linear dependence in speciation rate (anagenesis and cladogenesis) -# . ddmodel == 11 : linear dependence in speciation rate and immigration rate -# . ddmodel == 3 : linear dependence in extinction rate -# - cond = conditioning; currently only cond = 0 is possible -# . cond == 0 : no conditioning -# . cond == 1 : conditioning on presence on the island - - stop( - "This functionality is still under development and is not available yet." - ) - - 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) - idpars <- sort(c(idparsopt, idparsfix)) - missnumspec <- unlist(lapply(datalist, function(list) {list$missing_species})) - if (sum(missnumspec) > (res - 1)) { - cat("The number of missing species is too large relative to the resolution of the ODE.\n") - return(out2err) - } - if ((prod(idpars == (1:10)) != 1) || (length(initparsopt) != length(idparsopt)) || (length(parsfix) != length(idparsfix))) { - cat("The parameters to be optimized and/or fixed are incoherent.\n") - return(out2err) - } - if (length(idparsopt) > 10) { - cat("The number of parameters to be optimized is too high.\n") - return(out2err) - } - namepars <- c("area_pars1", "area_pars2", "area_pars3", "area_pars4", "lambda_c0", "mu_1", "mu_2", "K0", "gamma0", "lambda_a") - if (length(namepars[idparsopt]) == 0) { optstr = "nothing" } else { optstr = namepars[idparsopt] } - cat("You are optimizing", optstr, "\n") - if (length(namepars[idparsfix]) == 0) { fixstr = "nothing" } else { fixstr = namepars[idparsfix] } - cat("You are fixing", fixstr, "\n") - cat("Calculating the likelihood for the initial parameters.", "\n") - utils::flush.console() - trparsopt <- initparsopt / (1 + initparsopt) - trparsopt[which(initparsopt == Inf)] <- 1 - trparsfix <- parsfix / (1 + parsfix) - trparsfix[which(parsfix == Inf)] <- 1 - island_ontogeny <- translate_island_ontogeny(island_ontogeny) - pars2 <- c(res, ddmodel, cond, verbose, island_ontogeny, eqmodel = NA, tol, maxiter) - optimpars <- c(tol, maxiter) - initloglik <- DAISIE_loglik_all_choosepar3(trparsopt = trparsopt, trparsfix = trparsfix, idparsopt = idparsopt, idparsfix = idparsfix, pars2 = pars2, datalist = datalist, methode = methode, CS_version = CS_version, abstolint = tolint[1], reltolint = tolint[2]) - cat("The loglikelihood for the initial parameter values is", initloglik, "\n") - if (initloglik == -Inf) { - cat("The initial parameter values have a likelihood that is equal to 0 or below machine precision. Try again with different initial values.\n") - return(out2err) - } - cat("Optimizing the likelihood - this may take a while.", "\n") - utils::flush.console() - out <- - DDD::optimizer( - optimmethod = optimmethod, - optimpars = optimpars, - fun = DAISIE_loglik_all_choosepar3, - trparsopt = trparsopt, - idparsopt = idparsopt, - trparsfix = trparsfix, - idparsfix = idparsfix, - pars2 = pars2, - datalist = datalist, - methode = methode, - CS_version = CS_version, - abstolint = tolint[1], - reltolint = tolint[2], - jitter = jitter, - num_cycles = num_cycles - ) - if (out$conv != 0) { - cat("Optimization has not converged. Try again with different initial values.\n") - out2 <- out2err - out2$conv <- out$conv - return(out2) - } - MLtrpars <- as.numeric(unlist(out$par)) - MLpars <- MLtrpars / (1 - MLtrpars) - ML <- as.numeric(unlist(out$fvalues)) - MLpars1 <- rep(0, 10) - MLpars1[idparsopt] <- MLpars - if (length(idparsfix) != 0) { MLpars1[idparsfix] = parsfix } - if (MLpars1[8] > 10^7){ MLpars1[8] = Inf } - out2 <- data.frame(area_pars1 = MLpars1[1], area_pars2 = MLpars1[2], area_pars3 = MLpars1[3], area_pars4 = MLpars1[4], lambda_c0 = MLpars1[5], mu1 = MLpars1[6], mu2 = MLpars1[7], K0 = MLpars1[8], gamma0 = MLpars1[9], lambda_a = MLpars1[10], loglik = ML, df = length(initparsopt), conv = unlist(out$conv)) - s1 <- sprintf("Maximum likelihood parameter estimates: area_pars1: %f, area_pars2: %f, area_pars3: %f, area_pars4: %f, lambda_c0: %f, mu1: %f, mu2: %f, K0: %f, gamma0: %f, lambda_a: %f", MLpars1[1], MLpars1[2], MLpars1[3], MLpars1[4], MLpars1[5], MLpars1[6], MLpars1[7], MLpars1[8], MLpars1[9], MLpars1[10]) - s2 <- sprintf("Maximum loglikelihood: %f", ML) - cat("\n", s1, "\n", s2, "\n") - return(invisible(out2)) -} diff --git a/R/DAISIE_ML_CS.R b/R/DAISIE_ML_CS.R index 8db84dbe..933896f5 100644 --- a/R/DAISIE_ML_CS.R +++ b/R/DAISIE_ML_CS.R @@ -187,7 +187,7 @@ DAISIE_ML_CS <- DAISIE_ML <- function( if (datatype == "single") { if (is.na(island_ontogeny)) { - if(CS_version[[1]] == 2) { + if(CS_version[[1]] %in% c(2,3)) { out <- DAISIE_ML4(datalist = datalist, initparsopt = initparsopt, idparsopt = idparsopt, @@ -231,26 +231,10 @@ DAISIE_ML_CS <- DAISIE_ML <- function( jitter = jitter, num_cycles = num_cycles) } - } else - { - out <- DAISIE_ML3(datalist = datalist, - initparsopt = initparsopt, - idparsopt = idparsopt, - parsfix = parsfix, - idparsfix = idparsfix, - res = res, - ddmodel = ddmodel, - cond = cond, - island_ontogeny = island_ontogeny, - tol = tol, - maxiter = maxiter, - methode = methode, - optimmethod = optimmethod, - CS_version = CS_version, - verbose = verbose, - tolint = tolint, - jitter = jitter, - num_cycles = num_cycles) + } else { + stop( + "Time dependent estimation not yet available. Development ongoing." + ) } } else { diff --git a/R/DAISIE_SR_loglik_CS.R b/R/DAISIE_SR_loglik_CS.R index b83a7c00..adc596c0 100644 --- a/R/DAISIE_SR_loglik_CS.R +++ b/R/DAISIE_SR_loglik_CS.R @@ -32,7 +32,6 @@ odeproc <- function( } else if (times[1] < tshift & times[2] > tshift) { probs <- DAISIE_integrate(probs, c(times[1], tshift), fun, pars1, rtol = rtol, atol = atol, method = method) - y <- DAISIE_integrate(probs, c(tshift, times[2]), fun, pars2, rtol = rtol, atol = atol, method = method) } return(y) @@ -44,18 +43,25 @@ divdepvecproc <- function( k1, ddep, times, - type + lac_or_gam ) { tshift <- -abs(pars[11]) - if (type == "col") { - a <- 3 - } else { - a <- 0 - } if (times < tshift) { - return(divdepvec(pars[1 + a], pars[3], lx, k1, ddep)) + return(divdepvec( + lac_or_gam = lac_or_gam, + pars1 = pars[1:5], + lx = lx, + k1 = k1, + ddep = ddep + )) } else { - return(divdepvec(pars[6 + a], pars[8], lx, k1, ddep)) + return(divdepvec( + lac_or_gam = lac_or_gam, + pars1 = pars[6:10], + lx = lx, + k1 = k1, + ddep = ddep + )) } } @@ -99,7 +105,7 @@ DAISIE_SR_loglik_CS_M1 <- DAISIE_SR_loglik <- function( loglik <- -Inf return(loglik) } - if (sum(brts == 0) == 0) { + if (!any(brts == 0)) { brts[length(brts) + 1] <- 0 } # for stac = 0 and stac = 1, brts will contain origin of island and 0; length = 2; no. species should be 0 @@ -174,7 +180,7 @@ DAISIE_SR_loglik_CS_M1 <- DAISIE_SR_loglik <- function( k1 <- 1 } if (stac == 2 || stac == 3 || stac == 4) { - gamvec <- divdepvecproc(pars1, lx, k1, ddep * (ddep == 11 | ddep == 21), brts[2], "col") + gamvec <- divdepvecproc(pars1, lx, k1, ddep * (ddep == 11 | ddep == 21), brts[2], "gam") 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 @@ -190,7 +196,7 @@ DAISIE_SR_loglik_CS_M1 <- DAISIE_SR_loglik <- function( S1 <- length(brts) - 1 startk <- 3 if (S1 >= startk) { - lacvec <- divdepvecproc(pars1, lx, k1, ddep, brts[3], "spec") + lacvec <- divdepvecproc(pars1, lx, k1, ddep, brts[3], "lac") if (stac == 2 || stac == 3) { probs[1:lx] <- lacvec[1:lx] * (probs[1:lx] + probs[(2 * lx + 1):(3 * lx)]) probs[(lx + 1):(2 * lx)] <- lacvec[2:(lx + 1)] * probs[(lx + 1):(2 * lx)] @@ -209,11 +215,11 @@ DAISIE_SR_loglik_CS_M1 <- DAISIE_SR_loglik <- function( } for (k in startk:S1) { k1 <- k - 1 - probs <- odeproc(probs, brts[k:(k + 1)], DAISIE_loglik_rhs, c(pars1, k1, ddep), rtol = reltolint, atol = abstolint, method = methode) + probs <- odeproc(probs = probs, times = brts[k:(k + 1)], fun = DAISIE_loglik_rhs, pars = c(pars1, k1, ddep), rtol = reltolint, atol = abstolint, method = methode) cp <- checkprobs2(lx, loglik, probs, verbose); loglik <- cp[[1]]; probs <- cp[[2]] if (k < S1) { # speciation event - lacvec <- divdepvecproc(pars1, lx, k1, ddep, brts[k + 1], "spec") + lacvec <- divdepvecproc(pars1, lx, k1, ddep, brts[k + 1], "lac") probs[1:(2 * lx)] <- c(lacvec[1:lx], lacvec[2:(lx + 1)]) * probs[1:(2 * lx)] } } diff --git a/R/DAISIE_format_IW.R b/R/DAISIE_format_IW.R index 3cb9b44b..f9509aa4 100644 --- a/R/DAISIE_format_IW.R +++ b/R/DAISIE_format_IW.R @@ -102,8 +102,10 @@ DAISIE_format_IW_trait <- function(island_replicates, Mtotal <- M1 + M2 stt_all = matrix(ncol = 7,nrow = sample_freq + 1) colnames(stt_all) = c("Time","nI","nA","nC","nI2","nA2","nC2") - stt_all[,"Time"] = rev(seq(from = 0,to = total_time,length.out = sample_freq + 1)) - stt_all[1,2:7] = c(0,0,0,0,0,0) + + stt_all[,"Time"] = rev(seq(from = 0,to = total_time, length.out = sample_freq + 1)) + stt_all[1, 2:7] = c(0, 0, 0, 0, 0, 0) + the_stt = the_island$stt_table diff --git a/R/DAISIE_loglik_CS.R b/R/DAISIE_loglik_CS.R index 25607c7f..05949ee2 100644 --- a/R/DAISIE_loglik_CS.R +++ b/R/DAISIE_loglik_CS.R @@ -2,7 +2,7 @@ #' #' Sets or retrieves the max. number of iterations used by the odeint solver. #' -#' @param max_steps \code{num_threads}: sets max. iterations to \code{max_steps}. \cr +#' @param max_steps \code{max_steps}: sets max. iterations to \code{max_steps}. \cr #' @return current max. iterations #' #' @export DAISIE_CS_max_steps @@ -11,6 +11,18 @@ DAISIE_CS_max_steps <- function(max_steps) { } +#` adams_bashforth and adams_bashforth_moulton integration control +#' +#' Sets or retrieves the factor to calculate the step-size used by the odeint::adams_bashforth[_moulton] solvers. +#' +#' @param factor sets step-size to \code{factor * (t1 - t0)}. \cr +#' @return current factor +#' +#' @export DAISIE_abm_factor +DAISIE_abm_factor <- function(factor) { + return(.Call("daisie_odeint_abm_factor", factor)) +} + DAISIE_loglik_rhs_precomp <- function(pars,lx) { lac = pars[1] @@ -21,9 +33,9 @@ DAISIE_loglik_rhs_precomp <- function(pars,lx) kk = pars[6] ddep = pars[7] - nn = -2:(lx+2*kk+1) + nn = -2:(lx + 2 * kk + 1) lnn = length(nn) - nn = pmax(rep(0,lnn),nn) + nn = pmax(rep(0, lnn), nn) if(ddep == 0) { @@ -46,9 +58,21 @@ DAISIE_loglik_rhs_precomp <- function(pars,lx) } else if(ddep == 11) { laavec = laa * rep(1,lnn) - lacvec = pmax(rep(0,lnn),lac * (1 - nn/K)) + #lacvec = pmax(rep(0,lnn),lac * (1 - nn/K)) + lacvec <- get_clado_rate_per_capita(lac = lac, + d = 0, + num_spec = nn, + K = K, + A = 1) muvec = mu * rep(1,lnn) - gamvec = pmax(rep(0,lnn),gam * (1 - nn/K)) + muvec <- rep(1,lnn) * get_ext_rate_per_capita(mu = mu, + x = 0) + #gamvec = pmax(rep(0,lnn),gam * (1 - nn/K)) + gamvec <- get_immig_rate_per_capita(gam = gam, + num_spec = nn, + K = K, + A = 1) + } else if(ddep == 21) { laavec = laa * rep(1,lnn) @@ -66,6 +90,7 @@ DAISIE_loglik_rhs_precomp <- function(pars,lx) } DAISIE_loglik_rhs <- function(t, x, parsvec) { + rhs <- 0 kk <- parsvec[length(parsvec)] lx <- (length(x) - 1)/2 lnn <- lx + 4 + 2 * kk @@ -95,27 +120,90 @@ DAISIE_loglik_rhs <- function(t, x, parsvec) { ix3 = nil2lx ix4 = nil2lx-2 - dx1 = laavec[il1 + 1] * xx2[ix1] + + 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] + + -(muvec[il3 + 1] + lacvec[il3 + 1]) * nn[in3 + 1] * xx2[ix3] + + -laavec[il3 + 1] * xx2[ix3] + dx3 <- 0 + + return(list(c(dx1,dx2,dx3))) +} + +DAISIE_loglik_rhs1 <- function(t, x, parsvec) { + rhs <- 1 + kk <- parsvec[length(parsvec)] + lx <- (length(x))/4 + lnn <- lx + 4 + 2 * kk + laavec <- parsvec[1:lnn] + lacvec <- parsvec[(lnn + 1):(2 * lnn)] + muvec <- parsvec[(2 * lnn + 1):(3 * lnn)] + gamvec <- parsvec[(3 * lnn + 1):(4 * lnn)] + nn <- parsvec[(4 * lnn + 1):(5 * lnn)] + + xx1 <- c(0,0,x[1:lx],0) + xx2 <- c(0,0,x[(lx + 1):(2 * lx)],0) + xx3 <- c(0,0,x[(2 * lx + 1):(3 * lx)],0) + xx4 <- c(0,0,x[(3 * lx + 1):(4 * lx)],0) + + nil2lx <- 3:(lx + 2) + + il1 <- nil2lx+kk-1 + il2 <- nil2lx+kk+1 + il3 <- nil2lx+kk + il4 <- nil2lx+kk-2 + + in1 <- nil2lx+2*kk-1 + in2 <- nil2lx+1 + in3 <- nil2lx+kk + in4 <- nil2lx-1 + + ix1 <- nil2lx-1 + ix2 <- nil2lx+1 + ix3 <- nil2lx + ix4 <- nil2lx-2 + + dx1 <- lacvec[il1] * nn[in1] * xx1[ix1] + + laavec[il1 + 1] * xx2[ix1] + + lacvec[il4 + 1] * xx2[ix4] + + muvec[il2] * nn[in2] * xx1[ix2] + + muvec[il3 + 1] * xx2[ix3] + + -(muvec[il3] + lacvec[il3]) * nn[in3] * xx1[ix3] + + -gamvec[il3] * xx1[ix3] - dx2 = gamvec[il3] * xx1[ix3] + + dx2 <- gamvec[il3] * xx1[ix3] + + gamvec[il3] * xx3[ix3] + + gamvec[il3 + 1] * xx4[ix3] + 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] - dx3 = -(laavec[il3[1]] + lacvec[il3[1]] + gamvec[il3[1]] + muvec[il3[1]]) * xx3 + dx3 <- lacvec[il1] * nn[in1] * xx3[ix1] + + laavec[il1 + 1] * xx4[ix1] + + lacvec[il4 + 1] * xx4[ix4] + + muvec[il2] * nn[in2] * xx3[ix2] + + muvec[il3 + 1] * xx4[ix3] + + -(lacvec[il3] + muvec[il3]) * nn[in3] * xx3[ix3] + + -gamvec[il3] * xx3[ix3] - return(list(c(dx1,dx2,dx3))) + dx4 <- lacvec[il1 + 1] * nn[in1] * xx4[ix1] + + muvec[il2 + 1] * nn[in2] * xx4[ix2] + + -(lacvec[il3 + 1] + muvec[il3 + 1]) * nn[in3 + 1] * xx4[ix3] + + -gamvec[il3 + 1] * xx4[ix3] + + return(list(c(dx1,dx2,dx3,dx4))) } + DAISIE_loglik_rhs2 <- function(t, x, parsvec) { + rhs <- 2 kk <- parsvec[length(parsvec)] lx <- (length(x))/3 lnn <- lx + 4 + 2 * kk @@ -198,6 +286,7 @@ DAISIE_loglik_rhs2 <- function(t, x, parsvec) { return(list(c(dx1,dx2,dx3))) } + checkprobs <- function(lv, loglik, probs, verbose) { probs <- probs * (probs > 0) if (is.na(sum(probs[1:lv])) || is.nan(sum(probs))) { @@ -214,7 +303,7 @@ checkprobs <- function(lv, loglik, probs, verbose) { return(list(loglik, probs)) } -checkprobs2 <- function(lx, loglik, probs, verbose) { +checkprobs2 <- function(lv, loglik, probs, verbose) { probs <- probs * (probs > 0) if (is.na(sum(probs)) || is.nan(sum(probs))) { loglik <- -Inf @@ -230,20 +319,68 @@ checkprobs2 <- function(lx, loglik, probs, verbose) { return(list(loglik, probs)) } -divdepvec <- function(lacgam, +divdepvec <- function(lac_or_gam, pars1, + t = NULL, lx, k1, - ddep, - island_ontogeny = NA) { + ddep) { + island_ontogeny <- pars1[15] if (!is.na(island_ontogeny)) { - lacgamK <- divdepvec_time(lacgam, pars1, lx, k1, ddep, island_ontogeny) - lacgam <- lacgamK[1] - K <- lacgamK[2] + + # ddep is NOT yet used in the time dependent case + + # lac0 <- parsvec[1] + # mu0 <- parsvec[2] + # K0 <- parsvec[3] + # gam0 <- parsvec[4] + # laa0 <- parsvec[5] + # d <- parsvec[6] + # x <- parsvec[7] + # area_pars <- parsvec[8:14] + # island_ontogeny <- parsvec[15] + # sea_level <- parsvec[16] + # total_time <- parsvec[17] + # peak <- parsvec[18] + area <- island_area_vector( + timeval = abs(t), + area_pars = pars1[8:14], + island_ontogeny = island_ontogeny, + sea_level = pars1[16], + total_time = pars1[17], + peak = pars1[18] + ) + if (lac_or_gam == "lac") { + divdepvector <- get_clado_rate_per_capita( + lac = pars1[1], + d = pars1[6], + num_spec = (0:lx) + k1, + K = pars1[3], + A = area + ) + } else if (lac_or_gam == "gam") { + divdepvector <- get_immig_rate_per_capita( + gam = pars1[4], + num_spec = (0:lx) + k1, + K = pars1[3], + A = area + ) + } } else { - K <- pars1[1] + lacgam <- ifelse(lac_or_gam == "lac", pars1[1], pars1[4]) + divdepvector <- divdepvec1( + lacgam = lacgam, + K = pars1[3], + lx = lx, + k1 = k1, + ddep = ddep + ) + #divdepvector <- get_immig_rate_per_capita(gam = lacgam, + # num_spec = k1 + (0:lx), + # K = pars1[3], + # A = 1) } - return(divdepvec1(lacgam, K, lx, k1, ddep)) + return(divdepvector) } divdepvec1 <- function(lacgam, K, lx, k1, ddep) { @@ -294,46 +431,17 @@ DAISIE_loglik_CS_M1 <- DAISIE_loglik <- function(pars1, pars2[4] = 0 } ddep <- pars2[2] - cond <- pars2[3] - # TODO: check if pars2[5] should be NA of if this never happens - # if (is.na(pars2[5])) { - # pars2[5] <- 0 - # } - island_ontogeny <- pars2[5] - if (is.na(island_ontogeny)) # This calls the old code that doesn't expect - # ontogeny - { - lac <- pars1[1] - mu <- pars1[2] - K <- pars1[3] - if (ddep == 0) { - K <- Inf - } - gam <- pars1[4] - laa <- pars1[5] - pars1_in_divdepvec_call <- K + K <- pars1[3] + if (!is.na(pars2[5])) { + K <- K * pars1[8] + } + if(length(pars1) == 6) { + probability_of_init_presence <- pars1[6] + pars1 <- pars1[-6] } else { - #pars1[1:4] = area_pars - #pars1[5] = lac0 - #pars1[6:7] = mupars - #pars1[8] = K0 - #pars1[9] = gam0 - #pars1[10] = laa - #pars1[11] = island_ontogeny - pars1[11] <- island_ontogeny - - if (pars1[11] == 0 && pars1[6] != pars1[7]) { - warning("mu_min and mu_max are not equal! Setting mu_max = mu_min") - pars1[7] <- pars1[6] - } - - lac <- as.numeric(pars1[5]) - K <- as.numeric(pars1[8]) - gam <- as.numeric(pars1[9]) - pars1_in_divdepvec_call <- pars1 + probability_of_init_presence <- 0 } - - brts = -sort(abs(as.numeric(brts)),decreasing = TRUE) + brts <- -sort(abs(as.numeric(brts)),decreasing = TRUE) if(length(brts) == 1 & sum(brts == 0) == 1) { stop('The branching times contain only a 0. This means the island emerged at the present which is not allowed.'); @@ -389,6 +497,7 @@ DAISIE_loglik_CS_M1 <- DAISIE_loglik <- function(pars1, } N <- length(brts) - 1 # exception for N = 1 in high lambda case + lac <- pars1[1] if(lac == Inf & missnumspec == 0 & length(pars1) == 5 & N > 1) { loglik <- DAISIE_loglik_high_lambda(pars1, -brts, stac) } else { @@ -405,128 +514,162 @@ DAISIE_loglik_CS_M1 <- DAISIE_loglik <- function(pars1, # in all cases we integrate from the origin of the island to the colonization event # (stac 2, 3, 4), the first branching point (stac = 6, 7), to the maximum colonization # time (stac = 1, 5, 8, 9) or to the present (stac = 0) - probs = rep(0,2 * lx + 1) - probs[1] = 1 - k1 = 0 + probs <- rep(0,2 * lx + 1) + probs[1] <- 1 - probability_of_init_presence #Q^k_n + probs[lx + 1] <- probability_of_init_presence #Q^{M,k}_n + k1 <- 0 probs = DAISIE_integrate(probs,brts[1:2],DAISIE_loglik_rhs,c(pars1,k1,ddep),rtol = reltolint,atol = abstolint,method = methode) - cp = checkprobs(lv = 2 * lx, loglik, probs, verbose); loglik = cp[[1]]; probs = cp[[2]] + cp = checkprobs2(lv = 2 * lx, loglik, probs, verbose); loglik = cp[[1]]; probs = cp[[2]] if(stac == 0) - # for stac = 0, the integration is from the origin of the island until the present - # and we evaluate the probability of no clade being present and no immigrant species, - # but there can be missing species { - loglik = loglik + log(probs[1 + missnumspec]) - } else { - if (stac == 1 || stac == 5 || stac == 8 || stac == 9) - # for stac = 1, the integration is from the maximum - # colonization time (usually the island age + tiny time unit) - # until the present, where we set all probabilities where the - # immigrant is already present to 0 and we evaluate the - # probability of the immigrant species being present, - # but there can be missing species - # for stac = 5, we do exactly the same, but we evaluate the probability - # of an endemic species being present alone. - # for stac = 8 and 9, integration is from the maximum colonization - # time until the minimum colonization time + # for stac = 0, the integration was from the origin of the island until + # the present so we can immediately evaluate the probability of no clade + # being present and no immigrant species. + loglik = loglik + log(probs[1]) + } else + { + if (stac %in% c(1, 5:9) ) { - probs[(lx + 1):(2 * lx)] = 0 - probs = DAISIE_integrate(probs,brts[2:3],DAISIE_loglik_rhs,c(pars1,k1,ddep),rtol = reltolint,atol = abstolint,method = methode) - cp = checkprobs(lv = 2 * lx, loglik, probs, verbose); loglik = cp[[1]]; probs = cp[[2]] - if (stac == 1 || stac == 5) { - loglik = loglik + log(probs[(stac == 1) * lx + (stac == 5) + 1 + missnumspec]) - } else # stac = 8 or 9 - { - probs[(2 * lx + 1):(3 * lx)] = probs[1:lx] - probs[1:(2 * lx)] = 0 - k1 = 1 - probs = DAISIE_integrate(probs,c(brts[3:4]),DAISIE_loglik_rhs2,c(pars1,k1,ddep),rtol = reltolint,atol = abstolint,method = methode) - cp = checkprobs2(lx,loglik,probs, verbose); loglik = cp[[1]]; probs = cp[[2]] - loglik = loglik + log(probs[(stac == 8) * lx + (stac == 9) + 1 + missnumspec]) - } - } else { # stac = 2, 3, 4, 6, 7 - # for stac = 2, 3, 4, integration is then from the colonization event until the first branching time (stac = 2 and 3) or the present (stac = 4). We add a set of equations for Q_M,n, the probability that the process is compatible with the data, and speciation has not happened; during this time immigration is not allowed because it would alter the colonization time. - # After speciation, colonization is allowed again (re-immigration) - # all probabilities of states with the immigrant present are set to zero and all probabilities of states with endemics present are transported to the state with the colonist present waiting for speciation to happen. We also multiply by the (possibly diversity-dependent) immigration rate - # for stac = 6 and 7, integration is from the maximum colonization time until the first branching time - if (stac == 6 || stac == 7) { - probs[(lx + 1):(2 * lx)] = 0 - probs = DAISIE_integrate(probs,brts[2:3],DAISIE_loglik_rhs,c(pars1,k1,ddep),rtol = reltolint,atol = abstolint,method = methode) - cp = checkprobs(lv = 2 * lx, loglik, probs, verbose); loglik = cp[[1]]; probs = cp[[2]] - k1 = 1 + # for stac = 1, we now integrate from the maximum colonization time + # (usually the island age + tiny time unit) until the present, where + # we set all probabilities where the immigrant is already present to 0 + # and we evaluate the probability of the immigrant species being + # present, but there can be missing species. + # for stac = 5, we do exactly the same, but we evaluate the + # probability of an endemic species being present alone. + # for stac = 6 and 7, integration is from the maximum colonization + # time until the first branching time. This is the same as we did for + # stac = 1, 5. + # for stac = 8 and 9, integration is from the maximum colonization + # time until the minimum colonization time. + # In all cases we are dealing with a maximum colonization time which + # means that any colonization that took place before this maximum + # colonization time (including presence in the non-oceanic scenario) + # does not count and should be followed by another colonization. + # To allow this we introduce a third set of equations for the + # probability that colonization might have happened before but + # recolonization has not taken place yet (Q_M,n). + epss <- 1.01E-5 #We're taking the risk + if (abs(brts[2] - brts[1]) >= epss) { + probs[(2 * lx + 1):(4 * lx)] <- probs[1:(2 * lx)] + probs[1:(2 * lx)] <- 0 + } else { + probs[(2 * lx + 1):(4 * lx)] <- 0 } - 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) - 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 - probs = DAISIE_integrate(probs,c(brts[2:3]),DAISIE_loglik_rhs2,c(pars1,k1,ddep),rtol = reltolint,atol = abstolint,method = methode) - cp = checkprobs2(lx,loglik,probs, verbose); loglik = cp[[1]]; probs = cp[[2]] + + probs <- DAISIE_integrate(probs,brts[2:3],DAISIE_loglik_rhs1,c(pars1,k1,ddep),rtol = reltolint,atol = abstolint,method = methode) + cp <- checkprobs2(lx, loglik, probs, verbose); loglik <- cp[[1]]; probs <- cp[[2]] + if (stac %in% c(1, 5)) + { + loglik <- loglik + log(probs[(stac == 1) * lx + (stac == 5) + 1 + missnumspec]) + } else if (stac %in% c(6, 7, 8, 9)) + { + probs2 <- rep(0, 3 * lx) + probs2[1:(lx - 1)] <- (1:(lx - 1)) * probs[2:lx] + probs2[(lx + 1):(2 * lx - 1)] <- (1:(lx - 1)) * probs[(lx + 2):(2 * lx)] + probs2[2 * lx + 1] <- probs[(lx + 1)] + probs2[(2 * lx + 2):(3 * lx)] <- 0 + probs <- probs2 + rm(probs2) + if (stac %in% c(8, 9)) + { + k1 <- 1 + probs = DAISIE_integrate(probs,c(brts[3:4]),DAISIE_loglik_rhs2,c(pars1,k1,ddep),rtol = reltolint,atol = abstolint,method = methode) + cp = checkprobs2(lx, loglik, probs, verbose); loglik = cp[[1]]; probs = cp[[2]] + loglik = loglik + log(probs[(stac == 8) * (2 * lx + 1) + (stac == 9) + missnumspec]) + } } + } else if (stac %in% c(2, 3, 4) ) + { + # for stac = 2, 3, 4, integration is then from the colonization + # event until the first branching time (stac = 2 and 3) or the present + # (stac = 4). We add a set of equations for Q_M,n, the probability + # that the process is compatible with the data, and speciation has not + # happened; during this time immigration is not allowed because it + # would alter the colonization time. + t <- brts[2] + gamvec = divdepvec( + lac_or_gam = "gam", + pars1 = pars1, + t = t, + lx = lx, + k1 = k1, + ddep = ddep * (ddep == 11 | ddep == 21) + ) + 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 + probs = DAISIE_integrate(probs,c(brts[2:3]),DAISIE_loglik_rhs2,c(pars1,k1,ddep),rtol = reltolint,atol = abstolint,method = methode) + cp = checkprobs2(lx,loglik,probs, verbose); loglik = cp[[1]]; probs = cp[[2]] if (stac == 4) # if stac = 4, we're done and we take an element from Q_M,n { loglik = loglik + log(probs[2 * lx + 1 + missnumspec]) - } else { - # for stac = 2 and 3, at the first branching point all probabilities of states Q_M,n are transferred to probabilities where only endemics are present. Then go through the branching points. - S1 = length(brts) - 1 - startk = 3 - if(S1 >= startk) + } + } + if (stac %in% c(2, 3, 6, 7) ) + { + # at the first branching point all probabilities of states Q_M,n are + # transferred to probabilities where only endemics are present. Then + # go through the branching points. + S1 = length(brts) - 1 + startk = 3 + if(S1 >= startk) + { + t <- brts[startk] + lacvec <- divdepvec( + lac_or_gam = "lac", + pars1 = pars1, + t = t, + lx = lx + stac %in% c(6,7), + k1 = k1, + ddep = ddep + ) + if(stac %in% c(2,3)) { - t <- brts[startk] - lacvec <- divdepvec( - lac, - c(pars1_in_divdepvec_call, t, 1), - lx, - k1, - ddep, - island_ontogeny - ) - if (stac == 2 || stac == 3) { - probs[1:lx] <- lacvec[1:lx] * - (probs[1:lx] + probs[(2 * lx + 1):(3 * lx)]) - probs[(lx + 1):(2 * lx)] <- lacvec[2:(lx + 1)] * - probs[(lx + 1):(2 * lx)] - probs <- probs[-c((2 * lx + 2):(3 * lx))] - probs[2 * lx + 1] <- 0 - } - if(stac == 6 || stac == 7) - { - probs2 <- rep(0,2 * lx + 1) - probs2[(1:(lx - 1))] <- probs[(2:lx)] + 1/(2:lx) * probs[(lx + 1):(2 * lx - 1)] - probs2[lx] <- 1/(lx + 1) * probs[2 * lx] - probs2[(lx + 1):(2 * lx - 1)] <- (1:(lx - 1))/(2:lx) * probs[(lx + 2):(2 * lx)] - probs = probs2 - rm(probs2) - probs[1:lx] <- lacvec[1:lx] * probs[1:lx] - probs[(lx + 1):(2 * lx)] <- lacvec[2:(lx + 1)] * probs[(lx + 1):(2 * lx)] - } - for(k in startk:S1) + probs[1:lx] <- lacvec[1:lx] * (probs[1:lx] + probs[(2 * lx + 1):(3 * lx)]) + probs[(lx + 1):(2 * lx)] <- lacvec[2:(lx + 1)] * probs[(lx + 1):(2 * lx)] + } else { # stac in c(6,7) + probs[1:(lx - 1)] <- lacvec[2:lx] * + ((1:(lx - 1)) * probs[2:lx] + probs[(2 * lx + 1):(3 * lx - 1)]) + probs[(lx + 1):(2 * lx - 1)] <- lacvec[3:(lx + 1)] * (1:(lx - 1)) * + probs[(lx + 2):(2 * lx)] + probs[lx] <- 0 + probs[2 * lx] <- 0 + } + probs <- probs[-c((2 * lx + 2):(3 * lx))] + probs[2 * lx + 1] <- 0 + # After speciation, colonization is allowed again (re-immigration) + # all probabilities of states with the immigrant present are set to + # zero and all probabilities of states with endemics present are + # transported to the state with the colonist present waiting for + # speciation to happen. We also multiply by the (possibly diversity- + # dependent) immigration rate. + for (k in startk:S1) + { + k1 <- k - 1 + probs <- DAISIE_integrate(probs,brts[k:(k+1)],DAISIE_loglik_rhs,c(pars1,k1,ddep),rtol = reltolint,atol = abstolint,method = methode) + cp <- checkprobs2(lx, loglik, probs, verbose); loglik = cp[[1]]; probs = cp[[2]] + if(k < S1) { - k1 <- k - 1 - probs <- DAISIE_integrate(probs,brts[k:(k+1)],DAISIE_loglik_rhs,c(pars1,k1,ddep),rtol = reltolint,atol = abstolint,method = methode) - cp <- checkprobs2(lx, loglik, probs, verbose); loglik = cp[[1]]; probs = cp[[2]] - if(k < S1) - { - # speciation event - t <- brts[k + 1] - lacvec <- divdepvec( - lac, - c(pars1_in_divdepvec_call, t, 1), - lx, - k1, - ddep, - island_ontogeny - ) - probs[1:(2 * lx)] <- c(lacvec[1:lx], lacvec[2:(lx + 1)]) * - probs[1:(2 * lx)] - } + # speciation event + t <- brts[k + 1] + lacvec <- divdepvec( + lac_or_gam = "lac", + pars1 = pars1, + t = t, + lx = lx, + k1 = k1, + ddep = ddep + ) + probs[1:(2 * lx)] <- c(lacvec[1:lx], lacvec[2:(lx + 1)]) * + probs[1:(2 * lx)] } } - # we evaluate the probability of the phylogeny with any missing species at the present without (stac = 2 or stac = 6) or with (stac = 3 or stac = 7) the immigrant species - loglik <- loglik + log(probs[(stac == 3 || stac == 7) * lx + 1 + missnumspec]) } + # we evaluate the probability of the phylogeny with any missing species at the present without (stac = 2 or stac = 6) or with (stac = 3 or stac = 7) the immigrant species + loglik <- loglik + log(probs[(stac %in% c(3, 7)) * lx + 1 + missnumspec]) } } } @@ -552,7 +695,7 @@ DAISIE_loglik_CS_M1 <- DAISIE_loglik <- function(pars1, utils::flush.console() } if (is.na(loglik)) { - cat("NA in loglik encountered. Changing to -Inf.") + message("NA in loglik encountered. Changing to -Inf.") loglik <- -Inf } loglik <- as.numeric(loglik) @@ -561,17 +704,17 @@ DAISIE_loglik_CS_M1 <- DAISIE_loglik <- function(pars1, } DAISIE_loglik_CS_choice <- function( - pars1, - pars2, - datalist = NULL, - brts, - stac, - missnumspec, - methode = "lsodes", - CS_version = 1, - abstolint = 1E-16, - reltolint = 1E-10, - verbose = FALSE + pars1, + pars2, + datalist = NULL, + brts, + stac, + missnumspec, + methode = "lsodes", + CS_version = 1, + abstolint = 1E-16, + reltolint = 1E-10, + verbose = FALSE ) { if (CS_version[[1]] == 1) { @@ -724,20 +867,46 @@ approximate_logp0 <- function(gamma, mu, t) #' @export DAISIE_loglik_CS #' @export DAISIE_loglik_all DAISIE_loglik_CS <- DAISIE_loglik_all <- function( - pars1, - pars2, - datalist, - methode = "lsodes", - CS_version = 1, - abstolint = 1E-16, - reltolint = 1E-10) { + pars1, + pars2, + datalist, + methode = "lsodes", + CS_version = 1, + abstolint = 1E-16, + reltolint = 1E-10) { + if (length(pars1) == 14) { + if (datalist[[1]]$island_age > pars1[11]) { + stop( + "The island age in the area parameters is inconsistent with the island + data." + ) + } + peak <- calc_peak( + total_time = datalist[[1]]$island_age, + area_pars = create_area_pars( + max_area = pars1[8], + current_area = pars1[9], + proportional_peak_t = pars1[10], + total_island_age = pars1[11], + sea_level_amplitude = pars1[12], + sea_level_frequency = pars1[13], + island_gradient_angle = pars1[14] + ) + ) + pars1 <- c( + pars1, + island_ontogeny = pars2[5], + sea_level = pars2[6], + datalist[[1]]$island_age, + peak + ) + } pars1 <- as.numeric(pars1) cond <- pars2[3] endpars1 <- 5 - if(length(pars1) == 5 | !is.na(pars2[5])) # Normal no ont case - { + if(length(pars1) %in% c(5,6) | !is.na(pars2[5])) { if(!is.na(pars2[5])) { endpars1 <- length(pars1) @@ -755,7 +924,9 @@ DAISIE_loglik_CS <- DAISIE_loglik_all <- function( ) if(logp0 >= 0 & pars1[2]/pars1[1] > 100) { - logp0 <- approximate_logp0(gamma = pars1[4], mu = pars1[2], t = datalist[[1]]$island_age) + logp0 <- approximate_logp0(gamma = pars1[4], + mu = pars1[2], + t = datalist[[1]]$island_age) } if(logp0 >= 0) { @@ -786,7 +957,8 @@ DAISIE_loglik_CS <- DAISIE_loglik_all <- function( which(unlist(datalist)[which(names(unlist(datalist)) == "type1or2")] == 2) ) numimm_type1 <- length(datalist) - 1 - numimm_type2 - if (is.na(pars1[11]) == FALSE) { + + if (is.na(pars1[11]) == FALSE && length(pars1) == 11) { if (pars1[11] < numimm_type2 / numimm | pars1[11] > (1 - numimm_type1 / numimm)) { return(-Inf) @@ -839,7 +1011,7 @@ DAISIE_loglik_CS <- DAISIE_loglik_all <- function( numimm = c(datalist[[1]]$not_present_type1 + numimm_type1,datalist[[1]]$not_present_type2 + numimm_type2), logp0 = c(logp0_type1,logp0_type2) ) } - loglik = loglik - logcond + loglik <- loglik - logcond if(length(datalist) > 1) { @@ -922,19 +1094,61 @@ DAISIE_integrate_const <- function(initprobs,tvec,rhs_func,pars,rtol,atol,method # # Use a regular expression to extract if the part that we are interested # in is present - function_as_text <- as.character(body(rhs_func)[3]) - do_fun_1 <- grepl(pattern = "lx <- \\(length\\(x\\) - 1\\)/2", x = function_as_text) - do_fun_2 <- grepl(pattern = "lx <- \\(length\\(x\\)\\)/3", x = function_as_text) - if (do_fun_1) + function_as_text <- as.character(body(rhs_func)[2]) + if (function_as_text == 'rhs <- 0') { lx <- (length(initprobs) - 1)/2 parsvec <- c(DAISIE_loglik_rhs_precomp(pars,lx)) - y <- DAISIE_ode_cs(initprobs,tvec,parsvec,atol,rtol,method,runmod = "daisie_runmod") - } else if (do_fun_2) + y <- DAISIE_ode_cs( + initprobs, + tvec, + parsvec, + atol, + rtol, + method, + runmod = "daisie_runmod" + ) + #y <- deSolve::ode( + # y = initprobs, + # times = tvec, + # func = DAISIE_loglik_rhs1, + # parms = parsvec, + # rtol = rtol, + # atol = atol, + # method = method + # )[2, -1] + } else if (function_as_text == 'rhs <- 1') + { + lx <- (length(initprobs))/4 + parsvec <- c(DAISIE_loglik_rhs_precomp(pars,lx)) + y <- DAISIE_ode_cs(initprobs, + tvec, + parsvec, + atol, + rtol, + method, + runmod = "daisie_runmod1") + + } else if (function_as_text == 'rhs <- 2') { lx <- (length(initprobs))/3 parsvec <- c(DAISIE_loglik_rhs_precomp(pars,lx)) - y <- DAISIE_ode_cs(initprobs,tvec,parsvec,atol,rtol,method,runmod = "daisie_runmod2") + y <- DAISIE_ode_cs(initprobs, + tvec, + parsvec, + atol, + rtol, + method, + runmod = "daisie_runmod2") + #y <- deSolve::ode( + # y = initprobs, + # times = tvec, + # func = DAISIE_loglik_rhs2, + # parms = parsvec, + # rtol = rtol, + # atol = atol, + # method = method + #)[2, -1] } else { stop( @@ -948,46 +1162,61 @@ DAISIE_integrate_const <- function(initprobs,tvec,rhs_func,pars,rtol,atol,method # replacement for DAISIE_ode_FORTRAN # returns striped deSolve result DAISIE_ode_cs <- function( - initprobs, - tvec, - parsvec, - atol, - rtol, - methode, - runmod = "daisie_runmod") { + initprobs, + tvec, + parsvec, + atol, + rtol, + methode, + runmod = "daisie_runmod") { N <- length(initprobs) kk <- parsvec[length(parsvec)] if (runmod == "daisie_runmod") { lx <- (N - 1) / 2 + rhs_func <- DAISIE_loglik_rhs + } else if (runmod == "daisie_runmod1") + { + lx <- N / 4 + rhs_func <- DAISIE_loglik_rhs1 } else if (runmod == "daisie_runmod2") { lx <- N / 3 + rhs_func <- DAISIE_loglik_rhs2 } if (startsWith(methode, "odeint")) { probs <- .Call("daisie_odeint_cs", runmod, initprobs, tvec, lx, kk, parsvec[-length(parsvec)], methode, atol, rtol) - } - else { + } else if (startsWith(methode, "deSolve_R::")) { + methode <- substring(methode,12) + y <- deSolve::ode(y = initprobs, + times = tvec, + func = rhs_func, + parms = parsvec, + method = methode)[,1:(N + 1)] + probs <- y[-1,-1] + } else { y <- deSolve::ode(y = initprobs, parms = c(lx + 0.,kk + 0.), rpar = parsvec[-length(parsvec)], times = tvec, func = runmod, initfunc = "daisie_initmod", ynames = c("SV"), dimens = N + 2, nout = 1, outnames = c("Sum"), dllname = "DAISIE",atol = atol, rtol = rtol, method = methode)[,1:(N + 1)] - probs = y[-1,-1] # strip 1st row and 1st column + probs <- y[-1,-1] # strip 1st row and 1st column } return(probs) } # left in to cover the case this function is called from outside DAISIE_loglik_CS.R DAISIE_ode_FORTRAN <- function( - initprobs, - tvec, - parsvec, - atol, - rtol, - methode, - runmod = "daisie_runmod") { + initprobs, + tvec, + parsvec, + atol, + rtol, + methode, + runmod = "daisie_runmod") { N <- length(initprobs) kk <- parsvec[length(parsvec)] if (runmod == "daisie_runmod") { lx <- (N - 1) / 2 + } else if (runmod == "daisie_runmod1") { + lx <- N / 4 } else if (runmod == "daisie_runmod2") { lx <- N / 3 } diff --git a/R/DAISIE_loglik_CS_time.R b/R/DAISIE_loglik_CS_time.R index c4a652f2..0150f2a6 100644 --- a/R/DAISIE_loglik_CS_time.R +++ b/R/DAISIE_loglik_CS_time.R @@ -17,7 +17,12 @@ #' @references Valente, Luis M., Rampal S. Etienne, and Albert B. Phillimore. #' "The effects of island ontogeny on species diversity and phylogeny." #' Proceedings of the Royal Society of London B: Biological Sciences 281.1784 (2014): 20133227. -island_area_vector <- function(timeval, area_pars, island_ontogeny, sea_level) { +island_area_vector <- function(timeval, + area_pars, + island_ontogeny, + sea_level, + total_time, + peak) { # Constant if (island_ontogeny == 0 || is.na(island_ontogeny)) { if (area_pars[1] != 1 || is.null(area_pars[1])) { @@ -25,76 +30,76 @@ island_area_vector <- function(timeval, area_pars, island_ontogeny, sea_level) { } return(1) } else { # Ontogeny - area_pars <- create_area_pars(area_pars[1], - area_pars[2], - area_pars[3], - area_pars[4]) + area_pars <- create_area_pars(max_area = area_pars[1], + current_area = area_pars[2], + proportional_peak_t = area_pars[3], + total_island_age = area_pars[4], + sea_level_amplitude = area_pars[5], + sea_level_frequency = area_pars[6], + island_gradient_angle = area_pars[7]) area <- island_area( timeval = timeval, area_pars = area_pars, island_ontogeny = island_ontogeny, - sea_level = sea_level + sea_level = sea_level, + total_time = total_time, + peak = peak ) return(area) } } - - -#parsvec[1:4] = area_pars -#parsvec[5] = lac0 -#parsvec[6:7] = mupars -#parsvec[8] = K0 -#parsvec[9] = gam0 -#parsvec[10] = laa0 -#parsvec[11] = island_ontogeny -#parsvec[12] = kk -#parsvec[13] = ddep - DAISIE_loglik_rhs_time <- function(t, x, parsvec) { - kk <- parsvec[length(parsvec) - 1] + lac0 <- parsvec[1] + mu0 <- parsvec[2] + K0 <- parsvec[3] + gam0 <- parsvec[4] + laa0 <- parsvec[5] + d <- parsvec[6] + x_hyperpar <- parsvec[7] + area_pars <- parsvec[8:14] + island_ontogeny <- parsvec[15] + sea_level <- parsvec[16] + total_time <- parsvec[17] + peak <- parsvec[18] + kk <- parsvec[19] + ddep <- parsvec[20] + lx <- (length(x) - 1)/2 lnn <- lx + 4 + 2 * kk nn <- -2:(lx + 2 * kk + 1) - nn <- pmax(rep(0, lnn), nn) # Added this - area_pars <- parsvec[1:4] # to change - A <- parsvec[1] # to change - lac0 <- parsvec[5] - mu <- parsvec[6] - K0 <- parsvec[7] - gam0 <- parsvec[8] - laa0 <- parsvec[9] - island_ontogeny <- parsvec[10] - kk <- parsvec[11] - ddep <- parsvec[12] - time_for_area_calc <- abs(t) + nn <- pmax(rep(0, lnn), nn) + area <- island_area_vector( - timeval = time_for_area_calc, + timeval = abs(t), area_pars = area_pars, island_ontogeny = island_ontogeny, - sea_level = 0 + sea_level = sea_level, + total_time = total_time, + peak = peak ) - lac <- get_clado_rate( + lacvec <- get_clado_rate_per_capita( lac = lac0, - hyper_pars = create_hyper_pars(d = 0, - x = 0), - A = 0, # to change + d = d, + A = area, K = K0, - num_spec = 1 # Also need per capita?? + num_spec = nn + ) + muvec <- rep(1, lnn) * get_ext_rate_per_capita( + mu = mu0, + x = x_hyperpar, + A = area, + extcutoff = 1000000 ) - lacvec <- pmax(rep(0, lnn), lac0 * (1 - nn / (area * K0))) - mu <- get_ext_rate( - mu = mu, - hyper_pars = create_hyper_pars(d = 0, - x = 0), - A = A, # to change - extcutoff = 1100, - num_spec = 1 # Here we need per capita mu + gamvec <- get_immig_rate_per_capita( + gam = gam0, + A = area, + num_spec = nn, + K = K0 ) - muvec <- mu * rep(1, lnn) - gamvec <- pmax(rep(0, lnn), parsvec[9] * (1 - nn / (area * parsvec[8]))) - laavec <- parsvec[10] * rep(1, lnn) + laavec <- laa0 * rep(1, lnn) + xx1 <- c(0, 0, x[1:lx], 0) xx2 <- c(0, 0, x[(lx + 1):(2 * lx)], 0) xx3 <- x[2 * lx + 1] @@ -110,55 +115,83 @@ DAISIE_loglik_rhs_time <- function(t, x, parsvec) { ix2 <- nil2lx + 1 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] + -(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] + + + 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] + # The next two lines are relicts because the k = 1 case is dealth with by rhs2 + # 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] + -(muvec[il3 + 1] + lacvec[il3 + 1]) * nn[in3 + 1] * xx2[ix3] + -laavec[il3 + 1] * xx2[ix3] - dx3 <- -(laavec[il3[1]] + lacvec[il3[1]] + gamvec[il3[1]] + - muvec[il3[1]]) * xx3 + + # The next line is not relevant as xx3 is always 0 + # dx3 = -(laavec[il3[1]] + lacvec[il3[1]] + gamvec[il3[1]] + muvec[il3[1]]) * xx3 + # Still need to specify dx3 + dx3 <- 0 + return(list(c(dx1, dx2, dx3))) } -DAISIE_loglik_rhs_time2 <- function(t, x, parsvec) { - kk <- parsvec[length(parsvec) - 1] +DAISIE_loglik_rhs_time1 <- function(t, x, parsvec) { + lac0 <- parsvec[1] + mu0 <- parsvec[2] + K0 <- parsvec[3] + gam0 <- parsvec[4] + laa0 <- parsvec[5] + d <- parsvec[6] + x_hyperpar <- parsvec[7] + area_pars <- parsvec[8:14] + island_ontogeny <- parsvec[15] + sea_level <- parsvec[16] + total_time <- parsvec[17] + peak <- parsvec[18] + kk <- parsvec[19] + ddep <- parsvec[20] + lx <- (length(x))/3 lnn <- lx + 4 + 2 * kk nn <- -2:(lx + 2 * kk + 1) nn <- pmax(rep(0, lnn), nn) - area_pars <- parsvec[1:4] # to change - A <- parsvec[1] # to change - lac0 <- parsvec[5] - mu <- parsvec[6] - K0 <- parsvec[7] - gam0 <- parsvec[8] - laa0 <- parsvec[9] - island_ontogeny <- parsvec[10] - kk <- parsvec[11] - ddep <- parsvec[12] - time_for_area_calc <- abs(t) area <- island_area_vector( - timeval = time_for_area_calc, + timeval = abs(t), area_pars = area_pars, - island_ontogeny = island_ontogeny + island_ontogeny = island_ontogeny, + sea_level = sea_level, + total_time = total_time, + peak = peak ) - lacvec <- pmax(rep(0, lnn), lac0 * (1 - nn / (area * K0))) - mu <- get_ext_rate( - mu = mu, - hyper_pars = create_hyper_pars(d = 0, - x = 0), - A = A, - extcutoff = 1000, - num_spec = 1, # Here we need per capita mu + + lacvec <- get_clado_rate_per_capita( + lac = lac0, + d = d, + A = area, + K = K0, + num_spec = nn ) - muvec <- mu * rep(1, lnn) - gamvec <- pmax(rep(0, lnn), gam0 * (1 - nn / (area * K0))) + + muvec <- rep(1, lnn) * get_ext_rate_per_capita( + mu = mu0, + x = x_hyperpar, + A = area, + extcutoff = 1000000 + ) + gamvec <- get_immig_rate_per_capita( + gam = gam0, + A = area, + num_spec = nn, + K = K0 + ) + laavec <- laa0 * rep(1, lnn) xx1 <- c(0, 0, x[1:lx], 0) @@ -183,80 +216,181 @@ DAISIE_loglik_rhs_time2 <- function(t, x, parsvec) { ix4 <- nil2lx - 2 # inflow: - # anagenesis in colonist when k = 1: Q_M,n -> Q^1_n; n+k species present - # cladogenesis in colonist when k = 1: Q_M,n-1 -> Q^1_n; n+k-1 species - # present; rate twice + # recolonization when k = 0: Q_M,n -> Q^M,0_n + # rhs1 only applies to cases where k = 0 so this is actually not a relevant + # addition, but this indicates where rhs1 is critically different from rhs2. # anagenesis of reimmigrant: Q^M,k_n-1 -> Q^k,n; n+k-1+1 species present - # cladogenesis of reimmigrant: Q^M,k_n-2 -> Q^k,n; n+k-2+1 species present; - # rate once + # cladogenesis of reimmigrant: Q^M,k_n-2 -> Q^k,n; # extinction of reimmigrant: Q^M,k_n -> Q^k,n; n+k+1 species present # cladogenesis in one of the n+k-1 species: Q^k_n-1 -> Q^k_n; # n+k-1 species present; rate twice for k species - # extinction in one of the n+1 species: Q^k_n+1 -> Q^k_n; - # n+k+1 species present + # extinction in one of the n+1 species: Q^k_n+1 -> Q^k_n; n+k+1 species + # present # 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] + + 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] + -(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; n+k species present - # cladogenesis in n+k-1 species: Q^M,k_n-1 -> Q^M,k_n; n+k-1+1 species - # present; rate twice for k species + # immigration when there are n species: Q^0_M,n -> Q^M,0_n + # (This is where rhs1 is critically different from rhs2) + # immigration when there are n+k species: Q^k,n -> Q^M,k_n; + # n+k species present + # cladogenesis in n+k-1 species: Q^M,k_n-1 -> Q^M,k_n; + # n+k-1+1 species present; rate twice for k species # extinction in n+1 species: Q^M,k_n+1 -> Q^M,k_n; n+k+1+1 species present # 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] + + dx2 <- gamvec[il2 + 1] * xx3[ix3] * (kk == 0) + + gamvec[il2 + 1] * xx1[ix3] + + 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] - # only when k = 1 + # inflow: - # cladogenesis in one of the n-1 species: Q_M,n-1 -> Q_M,n; n+k-1 species present; rate once - # extinction in one of the n+1 species: Q_M,n+1 -> Q_M,n; n+k+1 species present + # cladogenesis in one of the n-1 species: Q_M,n-1 -> Q_M,n; + # n+k-1 species present; rate once + # extinction in one of the n+1 species: Q_M,n+1 -> Q_M,n; + # 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] + -(lacvec[il3] + muvec[il3]) * nn[in3] * xx3[ix3] + -(laavec[il3] + gamvec[il3]) * xx3[ix3] + return(list(c(dx1, dx2, dx3))) } -divdepvec_time <- function(lacgam, pars1, lx, k1, ddep, island_ontogeny) { - # pars1[1:4] = area_pars - # pars1[5] = lac0 - # pars1[6] = mu - # pars1[7] = K0 - # pars1[8] = gam0 - # pars1[9] = laa0 - # pars1[10] = island_ontogeny - # pars1[11] = t - # pars1[12] = 0 (for gam) or 1 (for lac) - area_pars <- pars1[1:4] - lac0 <- pars1[5] - mu <- pars1[6] - K0 <- pars1[7] - gam0 <- pars1[8] - laa0 <- pars1[9] - island_ontogeny <- pars1[10] - timeval <- pars1[11] - gamlac <- pars1[12] +DAISIE_loglik_rhs_time2 <- function(t, x, parsvec) { + lac0 <- parsvec[1] + mu0 <- parsvec[2] + K0 <- parsvec[3] + gam0 <- parsvec[4] + laa0 <- parsvec[5] + d <- parsvec[6] + x_hyperpar <- parsvec[7] + area_pars <- parsvec[8:14] + island_ontogeny <- parsvec[15] + sea_level <- parsvec[16] + total_time <- parsvec[17] + peak <- parsvec[18] + kk <- parsvec[19] + ddep <- parsvec[20] + + lx <- (length(x))/3 + lnn <- lx + 4 + 2 * kk + nn <- -2:(lx + 2 * kk + 1) + nn <- pmax(rep(0, lnn), nn) area <- island_area_vector( - timeval = timeval, + timeval = abs(t), area_pars = area_pars, - island_ontogeny = island_ontogeny + island_ontogeny = island_ontogeny, + sea_level = sea_level, + total_time = total_time, + peak = peak ) - lacA <- lac0 - gamA <- gam0 - KA <- area * K0 - lacgam <- (1 - gamlac) * gamA + gamlac * lacA - K <- KA - lacgamK <- c(lacgam, K) - return(lacgamK) + + lacvec <- get_clado_rate_per_capita( + lac = lac0, + d = d, + A = area, + K = K0, + num_spec = nn + ) + + muvec <- rep(1, lnn) * get_ext_rate_per_capita( + mu = mu0, + x = x_hyperpar, + A = area, + extcutoff = 1000000 + ) + gamvec <- get_immig_rate_per_capita( + gam = gam0, + A = area, + num_spec = nn, + K = K0 + ) + + laavec <- laa0 * rep(1, lnn) + + xx1 <- c(0, 0, x[1:lx], 0) + xx2 <- c(0, 0, x[(lx + 1):(2 * lx)], 0) + xx3 <- c(0, 0, x[(2 * lx + 1):(3 * lx)], 0) + + nil2lx <- 3:(lx + 2) + + il1 <- nil2lx + kk - 1 + il2 <- nil2lx + kk + 1 + il3 <- nil2lx + kk + il4 <- nil2lx + kk - 2 + + in1 <- nil2lx + 2 * kk - 1 + in2 <- nil2lx + 1 + in3 <- nil2lx + kk + in4 <- nil2lx - 1 + + ix1 <- nil2lx - 1 + ix2 <- nil2lx + 1 + ix3 <- nil2lx + ix4 <- nil2lx - 2 + + # inflow: + # anagenesis in colonist when k = 1: Q_M,n -> Q^1_n; n+k species present + # cladogenesis in colonist when k = 1: Q_M,n-1 -> Q^1_n; + # n+k-1 species present; rate twice + # anagenesis of reimmigrant: Q^M,k_n-1 -> Q^k,n; n+k-1+1 species present + # cladogenesis of reimmigrant: Q^M,k_n-2 -> Q^k,n; + # n+k-2+1 species present; rate once + # extinction of reimmigrant: Q^M,k_n -> Q^k,n; n+k+1 species present + # cladogenesis in one of the n+k-1 species: Q^k_n-1 -> Q^k_n; + # n+k-1 species present; rate twice for k species + # extinction in one of the n+1 species: Q^k_n+1 -> Q^k_n; n+k+1 species + # present + # 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] + + # inflow: + # immigration when there are n+k species: Q^k,n -> Q^M,k_n; + # n+k species present + # cladogenesis in n+k-1 species: Q^M,k_n-1 -> Q^M,k_n; + # n+k-1+1 species present; rate twice for k species + # extinction in n+1 species: Q^M,k_n+1 -> Q^M,k_n; n+k+1+1 species present + # 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] + + -(muvec[il3 + 1] + lacvec[il3 + 1]) * nn[in3 + 1] * xx2[ix3] + + -laavec[il3 + 1] * xx2[ix3] + + # only when k = 1 + # inflow: + # cladogenesis in one of the n-1 species: Q_M,n-1 -> Q_M,n; + # n+k-1 species present; rate once + # extinction in one of the n+1 species: Q_M,n+1 -> Q_M,n; + # 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] + + -(lacvec[il3] + muvec[il3]) * nn[in3] * xx3[ix3] + + -(laavec[il3] + gamvec[il3]) * xx3[ix3] + + return(list(c(dx1, dx2, dx3))) } DAISIE_integrate_time <- function(initprobs, @@ -266,10 +400,11 @@ DAISIE_integrate_time <- function(initprobs, rtol, atol, method) { - function_as_text <- as.character(body(rhs_func)[3]) - do_fun_1 <- grepl(pattern = "lx <- \\(length\\(x\\) - 1\\)/2", x = function_as_text) - do_fun_2 <- grepl(pattern = "lx <- \\(length\\(x\\)\\)/3", x = function_as_text) - if (do_fun_1) { + function_as_text <- as.character(body(rhs_func)[2]) + do_fun <- grepl(pattern = "rhs <- 0",x = function_as_text) + do_fun_1 <- grepl(pattern = "rhs <- 1",x = function_as_text) + do_fun_2 <- grepl(pattern = "rhs <- 2",x = function_as_text) + if (do_fun) { y <- deSolve::ode( initprobs, tvec, @@ -279,6 +414,16 @@ DAISIE_integrate_time <- function(initprobs, rtol = rtol, method = method ) + } else if (do_fun_1) { + y <- deSolve::ode( + initprobs, + tvec, + func = DAISIE_loglik_rhs_time1, + pars, + atol = atol, + rtol = rtol, + method = method + ) } else if (do_fun_2) { y <- deSolve::ode( initprobs, @@ -295,5 +440,6 @@ DAISIE_integrate_time <- function(initprobs, "Value of 'function_as_text':", function_as_text ) } + y <- y[-1,-1] return(y) } diff --git a/R/DAISIE_loglik_IW.R b/R/DAISIE_loglik_IW.R index f09083fd..e55b5c94 100644 --- a/R/DAISIE_loglik_IW.R +++ b/R/DAISIE_loglik_IW.R @@ -259,6 +259,8 @@ DAISIE_loglik_rhs_IW <- function(t,x,cp) #' \code{'odeint::runge_kutta_fehlberg78'} [default] #' \code{'odeint::runge_kutta_dopri5'} #' \code{'odeint::bulirsch_stoer'} +#' \code{'odeint::adams_bashforth_[1|2|3|4|5|6|7|8]} +#' \code{'odeint::adams_bashforth_moulton_[1|2|3|4|5|6|7|8]} #' without \code{odeint::}-prefix, \code{\link[deSolve]{ode}} method is #' assumed. #' @param abstolint Absolute tolerance of the integration diff --git a/R/DAISIE_loglik_high_lambda.R b/R/DAISIE_loglik_high_lambda.R index 6673ddea..ccbb2d8b 100644 --- a/R/DAISIE_loglik_high_lambda.R +++ b/R/DAISIE_loglik_high_lambda.R @@ -24,4 +24,4 @@ DAISIE_loglik_high_lambda <- function(pars1, brts, stac) { out <- -Inf } return(out) -} \ No newline at end of file +} diff --git a/R/DAISIE_max_rates.R b/R/DAISIE_max_rates.R index 4e01285f..45e6857e 100644 --- a/R/DAISIE_max_rates.R +++ b/R/DAISIE_max_rates.R @@ -60,6 +60,8 @@ update_max_rates <- function(gam, ana_max_rate = ana_max_rate, clado_max_rate = clado_max_rate ) + # print("max rates") + # print(max_rates) return(max_rates) } diff --git a/R/DAISIE_rates.R b/R/DAISIE_rates.R index 45a235de..4f44da80 100644 --- a/R/DAISIE_rates.R +++ b/R/DAISIE_rates.R @@ -83,6 +83,7 @@ update_rates <- function(timeval, K = K, mainland_n = mainland_n ) + # testit::assert(is.numeric(immig_rate)) ext_rate <- get_ext_rate( mu = mu, @@ -91,6 +92,7 @@ update_rates <- function(timeval, num_spec = num_spec, A = A ) + # testit::assert(is.numeric(ext_rate)) ana_rate <- get_ana_rate( laa = laa, @@ -104,6 +106,7 @@ update_rates <- function(timeval, K = K, A = A ) + # testit::assert(is.numeric(clado_rate)) rates <- list( @@ -278,10 +281,41 @@ island_area <- function(timeval, } } -#' Function to describe changes in extinction rate through time. +#' Function to describe per-capita changes in extinction rate through time +#' +#' This function is only called directly inside the RHS of the ontogeny +#' likelihood functions. In all other cases \code{\link{get_ext_rate}()} is to +#' be called instead. #' #' @inheritParams default_params_doc #' +#' @return Numeric with per capita extinction rate, given A(t), x, and mu0. +#' @keywords internal +#' +#' @examples +#' ext_rate_per_capita <- DAISIE:::get_ext_rate_per_capita( +#' mu = 0.5, +#' x = 1, +#' A = 1000 +#' ) +get_ext_rate_per_capita <- function(mu, + x, + extcutoff = 1000, + A = 1) { + ext_rate_per_capita <- max(0, mu * (A ^ -x), na.rm = TRUE) + ext_rate_per_capita <- min(ext_rate_per_capita, extcutoff, na.rm = TRUE) + return(ext_rate_per_capita) +} + + +#' Calculate extinction rate +#' +#' @inheritParams default_params_doc +#' +#' @return A numeric, with the extinction rate given the base extinction rate, +#' if present, the hyperparemeter \code{x}, A(t) if time-dependent and traits +#' if running a traint model. +#' #' @keywords internal #' @family rate calculations #' @references Valente, Luis M., Rampal S. Etienne, and Albert B. Phillimore. @@ -293,13 +327,19 @@ get_ext_rate <- function(mu, hyper_pars, extcutoff = 1000, num_spec, - A, + A = 1, trait_pars = NULL, island_spec = NULL) { x <- hyper_pars$x if (is.null(trait_pars)) { - ext_rate <- max(0, mu * (A ^ -x) * num_spec, na.rm = TRUE) + + ext_rate <- num_spec * get_ext_rate_per_capita( + mu = mu, + x = x, + extcutoff = extcutoff, + A = A + ) ext_rate <- min(ext_rate, extcutoff, na.rm = TRUE) # testit::assert(ext_rate >= 0) return(ext_rate) @@ -320,6 +360,8 @@ get_ext_rate <- function(mu, } } + + #' Calculate anagenesis rate #' @description Internal function. #' Calculates the anagenesis rate given the current number of @@ -341,11 +383,13 @@ get_ana_rate <- function(laa, # testit::assert(is.numeric(ana_rate)) # testit::assert(ana_rate >= 0) return(ana_rate) - }else{ + } else { ana_rate1 = laa * length(intersect(which(island_spec[,4] == "I"), which(island_spec[,8] == "1"))) - ana_rate2 = trait_pars$ana_rate2 * length(intersect(which(island_spec[,4] == "I"), - which(island_spec[,8] == "2"))) + ana_rate2 = trait_pars$ana_rate2 * length( + intersect(which(island_spec[,4] == "I"), + which(island_spec[,8] == "2")) + ) # testit::assert(is.numeric(ana_rate1)) # testit::assert(ana_rate1 >= 0) # testit::assert(is.numeric(ana_rate2)) @@ -356,6 +400,36 @@ get_ana_rate <- function(laa, } } +#' Calculate per-capita cladogenesis rate +#' +#' This function is only called directly inside the RHS of the ontogeny +#' likelihood functions. In all other cases \code{\link{get_clado_rate}()} is to +#' be called instead. +#' +#' @inheritParams default_params_doc +#' +#' @return Numeric with the per-capita cladogenesis rate given a base +#' cladogenesis rate, K, A and the d hyperparameter. +#' @keywords internal +#' +#' @examples +#' lac <- 0.4 +#' d <- 0 +#' num_spec <- 2 +#' K <- 10 +#' A <- 1 +#' clado_rate_pc <- DAISIE:::get_clado_rate_per_capita(lac, d, num_spec, K, A) +get_clado_rate_per_capita <- function(lac, + d, + num_spec, + K, + A = 1) { + clado_rate_per_capita <- lac * (A ^ d) * (1 - num_spec / (K * A)) + clado_rate_per_capita <- pmax(0, clado_rate_per_capita, na.rm = TRUE) + + return(clado_rate_per_capita) +} + #' Calculate cladogenesis rate #' @description Internal function. #' Calculates the cladogenesis rate given the current number of @@ -377,22 +451,39 @@ get_clado_rate <- function(lac, d <- hyper_pars$d if (is.null(trait_pars)) { - clado_rate <- max( - 0, lac * num_spec * (A ^ d) * (1 - num_spec / (K * A)), na.rm = TRUE - ) - # testit::assert(clado_rate >= 0) - # testit::assert(is.numeric(clado_rate)) - return(clado_rate) - }else{ + clado_rate_pc <- get_clado_rate_per_capita( + lac = lac, + d = d, + num_spec = num_spec, + K = K, + A = A + ) + clado_rate <- num_spec * clado_rate_pc + # testit::assert(clado_rate >= 0) + # testit::assert(is.numeric(clado_rate)) + return(clado_rate) + } else { num_spec_trait1 <- length(which(island_spec[, 8] == "1")) num_spec_trait2 <- length(which(island_spec[, 8] == "2")) - clado_rate1 <- max( + clado_rate1_pc <- max( 0, lac * num_spec_trait1 * (1 - num_spec / K), na.rm = TRUE) - clado_rate2 <- max( - 0, trait_pars$clado_rate2 * num_spec_trait2 * (1 - num_spec / K), - na.rm = TRUE + clado_rate1_pc <- get_clado_rate_per_capita( + lac = lac, + d = d, + num_spec = num_spec, + K = K, + A = A + ) + clado_rate1 <- num_spec_trait1 * clado_rate1_pc + clado_rate2_pc <- get_clado_rate_per_capita( + lac = lac, + d = d, + num_spec = num_spec_trait2, + K = K, + A = A ) + clado_rate2 <- num_spec_trait2* clado_rate2_pc # testit::assert(clado_rate1 >= 0) # testit::assert(clado_rate2 >= 0) # testit::assert(is.numeric(clado_rate1)) @@ -403,6 +494,34 @@ get_clado_rate <- function(lac, } } +#' Calculate per-capita immigration rate +#' +#' This function is only called directly inside the RHS of the ontogeny +#' likelihood functions. In all other cases \code{\link{get_immig_rate}()} is to +#' be called instead. +#' +#' @inheritParams default_params_doc +#' +#' @return A numeric with the per-capita immigration rate given A(t) and K. +#' @keywords internal +#' +#' @examples +#' immig_rate_per_capita <- DAISIE:::get_immig_rate_per_capita( +#' gam = 0.001, +#' num_spec = 5, +#' K = 20, +#' A = 1000 +#' ) +get_immig_rate_per_capita <- function(gam, + num_spec, + K, + A = 1) { + immig_rate_per_capita <- pmax( + 0, gam * (1 - (num_spec / (K * A))), na.rm = TRUE + ) + return(immig_rate_per_capita) +} + #' Calculate immigration rate #' @description Internal function. #' Calculates the immigration rate given the current number of @@ -417,7 +536,7 @@ get_clado_rate <- function(lac, #' "The effects of island ontogeny on species diversity and phylogeny." #' Proceedings of the Royal Society of London B: Biological Sciences 281.1784 (2014): 20133227. get_immig_rate <- function(gam, - A, + A = 1, num_spec, K, mainland_n, @@ -425,18 +544,30 @@ get_immig_rate <- function(gam, island_spec = NULL) { if (is.null(trait_pars)) { - immig_rate <- max(c(mainland_n * gam * (1 - (num_spec / (A * K))), - 0), na.rm = TRUE) + immig_rate <- mainland_n * get_immig_rate_per_capita( + gam = gam, + num_spec = num_spec, + K = K, + A = A + ) # testit::assert(is.numeric(immig_rate)) # testit::assert(immig_rate >= 0) return(immig_rate) } else { mainland_n2 <- trait_pars$M2 gam2 <- trait_pars$immig_rate2 - immig_rate1 <- max(c(mainland_n * gam * (1 - (num_spec / (A * K))), - 0), na.rm = TRUE) - immig_rate2 <- max(c(mainland_n2 * gam2 * (1 - (num_spec / (A * K))), - 0), na.rm = TRUE) + immig_rate1 <- max(0, mainland_n * get_immig_rate_per_capita( + gam = gam, + num_spec = num_spec, + K = K, + A = A + ), na.rm = TRUE) + immig_rate2 <- max(0, mainland_n2 * get_immig_rate_per_capita( + gam = gam2, + num_spec = num_spec, + K = K, + A = A + ), na.rm = TRUE) # testit::assert(is.numeric(immig_rate1)) # testit::assert(immig_rate1 >= 0) # testit::assert(is.numeric(immig_rate2)) @@ -469,20 +600,20 @@ get_immig_rate <- function(gam, get_trans_rate <- function(trait_pars, island_spec){ - # Make function accept island_spec matrix or numeric - if (is.matrix(island_spec) || is.null(island_spec)) { - num_spec_trait1 <- length(which(island_spec[, 8] == "1")) - num_spec_trait2 <- length(which(island_spec[, 8] == "2")) - } - trans_rate1 <- trait_pars$trans_rate * num_spec_trait1 - trans_rate2 <- trait_pars$trans_rate2 * num_spec_trait2 - testit::assert(is.numeric(trans_rate1)) - testit::assert(trans_rate1 >= 0) - testit::assert(is.numeric(trans_rate2)) - testit::assert(trans_rate2 >= 0) - trans_list <- list(trans_rate1 = trans_rate1, - trans_rate2 = trans_rate2) - return(trans_list) + # Make function accept island_spec matrix or numeric + if (is.matrix(island_spec) || is.null(island_spec)) { + num_spec_trait1 <- length(which(island_spec[, 8] == "1")) + num_spec_trait2 <- length(which(island_spec[, 8] == "2")) + } + trans_rate1 <- trait_pars$trans_rate * num_spec_trait1 + trans_rate2 <- trait_pars$trans_rate2 * num_spec_trait2 + testit::assert(is.numeric(trans_rate1)) + testit::assert(trans_rate1 >= 0) + testit::assert(is.numeric(trans_rate2)) + testit::assert(trans_rate2 >= 0) + trans_list <- list(trans_rate1 = trans_rate1, + trans_rate2 = trans_rate2) + return(trans_list) } @@ -502,9 +633,11 @@ calc_next_timeval <- function(max_rates, timeval) { # testit::assert(timeval >= 0) if (length(max_rates) == 4) { ## no considering about two trait states - totalrate <- max_rates[[1]] + max_rates[[2]] + max_rates[[3]] + max_rates[[4]] + totalrate <- + max_rates[[1]] + max_rates[[2]] + max_rates[[3]] + max_rates[[4]] } else { - totalrate <- max_rates[[1]] + max_rates[[2]] + max_rates[[3]] + max_rates[[4]] + + totalrate <- + max_rates[[1]] + max_rates[[2]] + max_rates[[3]] + max_rates[[4]] + max_rates[[5]] + max_rates[[6]] + max_rates[[7]] + max_rates[[8]] + max_rates[[9]] + max_rates[[10]] } @@ -567,7 +700,7 @@ calc_Abeta <- function(proptime, a <- f * peak / (1 + f) b <- peak / (1 + f) At <- Amax * proptime ^ a * - (1 - proptime) ^ b / ((a / (a + b)) ^ a * (b / (a + b)) ^ b) + (1 - proptime) ^ b / ((a / (a + b)) ^ a * (b / (a + b)) ^ b) return(At) } diff --git a/R/DAISIE_sim_core_time_dep.R b/R/DAISIE_sim_core_time_dep.R index 92e85d9a..a3c988db 100644 --- a/R/DAISIE_sim_core_time_dep.R +++ b/R/DAISIE_sim_core_time_dep.R @@ -67,7 +67,6 @@ DAISIE_sim_core_time_dep <- function( num_spec <- length(island_spec[, 1]) num_immigrants <- length(which(island_spec[, 4] == "I")) - #### Start Monte Carlo #### while (timeval < total_time) { max_rates <- update_max_rates( @@ -110,6 +109,10 @@ DAISIE_sim_core_time_dep <- function( num_immigrants = num_immigrants, mainland_n = mainland_n ) + # print("rates") + # print(rates) + # print(island_spec) + # print(timeval) # testit::assert(are_rates(rates)) possible_event <- DAISIE_sample_event_time_dep( max_rates = max_rates diff --git a/R/DAISIE_utils.R b/R/DAISIE_utils.R index c7a6138a..8c34aa30 100644 --- a/R/DAISIE_utils.R +++ b/R/DAISIE_utils.R @@ -432,3 +432,24 @@ DAISIE_spec_tables <- function(stt_table, island_spec = island_spec, maxspecID = maxspecID)) } + +#' Add a column to a data frame +#' +#' @param df data frame to add the column to +#' @param position location in data frame where to insert the column. +#' Position can also be a name of a column +#' @param column_to_insert the elements of the column to insert. If +#' the column has a name, this name will be copied into the data frame. +#' Id is does not have a name, it will get the name "nc". +#' +#' @return A data frame with the column inserted +add_column_to_dataframe <- function(df, position, column_to_insert) { + if(is.character(position)) { + position <- which(names(df) == position) + } + df <- data.frame(df[1:position], + nc = column_to_insert, + df[(position + 1):ncol(df)]) + names(df)[names(df) == 'nc'] <- names(column_to_insert) + return(df) +} diff --git a/R/default_params_doc.R b/R/default_params_doc.R index 334a6051..f8381be2 100644 --- a/R/default_params_doc.R +++ b/R/default_params_doc.R @@ -105,7 +105,9 @@ #' \code{\link{is_island_ontogeny_input}()}. \cr In all other functions a #' numeric describing the type of island ontogeny. Can be \code{0} for #' constant, \code{1} for a beta function describing area through time. In ML -#' functions \code{island_ontogeny = NA} assumes constant ontogeny. +#' functions \code{island_ontogeny = NA} assumes constant ontogeny. Time +#' dependent estimation is not yet available as development is still ongoing. +#' Will return an error if called in that case. #' @param sea_level In \code{\link{DAISIE_sim_time_dep}()} and plotting a #' string describing the type of sea level. Can be \code{"const"} or #' \code{"sine"} for a sine function describing area through time. String @@ -268,7 +270,7 @@ #' \code{"odeint::runge_kutta_fehlberg78"} #' \code{"odeint::runge_kutta_dopri5"} #' \code{"odeint::bulirsch_stoer"} -#' without \code{odeint::}-prefix, \code{\link[deSolve]{ode}} method is +#' without \code{odeint::}-prefix, \code{\link{deSolve}{ode}} method is #' assumed. The default method overall is #' \code{"lsodes"} for \code{\link{DAISIE_ML_CS}()} #' and \code{"ode45"} from \code{\link[deSolve]{ode}()} for @@ -413,7 +415,7 @@ #' at the end of the simulation). #' @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 +#' 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. diff --git a/README.md b/README.md index ab55205f..e01d37ef 100644 --- a/README.md +++ b/README.md @@ -24,7 +24,7 @@ The model can be fitted to both empirical dated phylogenies and simulated data. **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.** -The DAISIE package has a stable version on [CRAN](https://CRAN.R-project.org/package=DAISIE) and a development version on GitHub. +The `DAISIE` package has a stable version on [CRAN](https://CRAN.R-project.org/package=DAISIE) and a development version on GitHub. ### From CRAN @@ -37,7 +37,7 @@ install.packages("DAISIE") ### From GitHub -Install DAISIE from this GitHub repository by running: +Install `DAISIE` from this GitHub repository by running: ``` r install.packages("remotes") @@ -71,7 +71,11 @@ Remotes: rsetienne/DAISIE ``` -## `git` branching model +## Support/Questions + +For feature requests or bug-reports or other matters, please submit an [issue](https://github.com/rsetienne/DAISIE/issues/new). + +## `git` branching workflow * `master`: build should always pass. [@rsetienne](https://github.com/rsetienne) has control over `develop` to `master` merges. * `develop`: merge of topic branches, merge with `master` by [@rsetienne](https://github.com/rsetienne) iff build passes. @@ -100,9 +104,10 @@ Valente, L., Phillimore, A.B., Melo, M., Warren, B.H., Clegg, S.M., Havenstein, 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. EcoEvoRxiv. https://doi.org/10.32942/osf.io/w5ntf +Valente, L., Kristensen, N., Phillimore, A. B., & Etienne, R. S. (2021). Report of programming bugs in the DAISIE R package: consequences and correction. https://doi.org/10.32942/osf.io/w5ntf + -Santos Neves, P., Lambert, J. W., Valente, L., & Etienne, R. S. (2021).The robustness of a simple dynamic model of island biodiversity to geological and eustatic change. bioRxiv 2021.07.26.453064. https://doi.org/10.1101/2021.07.26.453064 +Santos Neves, P.\*, Lambert, J. W.\*, Valente, L., & Etienne, R. S. (2021). The robustness of a simple dynamic model of island biodiversity to geological and eustatic change. bioRxiv. https://doi.org/10.1101/2021.07.26.453064 Lambert, J. W., Santos Neves, P., Bilderbeek, R. L. C., Valente, L., Etienne, R. S. (2022). The effect of mainland dynamics on data and parameter estimates in island biogeography. bioRxiv. https://doi.org/10.1101/2022.01.13.476210 diff --git a/man/DAISIE_CS_max_steps.Rd b/man/DAISIE_CS_max_steps.Rd index 7419ce30..92a38c87 100644 --- a/man/DAISIE_CS_max_steps.Rd +++ b/man/DAISIE_CS_max_steps.Rd @@ -7,7 +7,7 @@ DAISIE_CS_max_steps(max_steps) } \arguments{ -\item{max_steps}{\code{num_threads}: sets max. iterations to \code{max_steps}. \cr} +\item{max_steps}{\code{max_steps}: sets max. iterations to \code{max_steps}. \cr} } \value{ current max. iterations diff --git a/man/DAISIE_ML.Rd b/man/DAISIE_ML.Rd index 7bba693c..057774ac 100644 --- a/man/DAISIE_ML.Rd +++ b/man/DAISIE_ML.Rd @@ -126,7 +126,9 @@ describing area through time. String checked by \code{\link{is_island_ontogeny_input}()}. \cr In all other functions a numeric describing the type of island ontogeny. Can be \code{0} for constant, \code{1} for a beta function describing area through time. In ML -functions \code{island_ontogeny = NA} assumes constant ontogeny.} +functions \code{island_ontogeny = NA} assumes constant ontogeny. Time +dependent estimation is not yet available as development is still ongoing. +Will return an error if called in that case.} \item{eqmodel}{Sets the equilibrium constraint that can be used during the likelihood optimization. Only available for datatype = 'single'.\cr\cr @@ -159,7 +161,7 @@ solvers (steppers) are: \code{"odeint::runge_kutta_fehlberg78"} \code{"odeint::runge_kutta_dopri5"} \code{"odeint::bulirsch_stoer"} -without \code{odeint::}-prefix, \code{\link[deSolve]{ode}} method is +without \code{odeint::}-prefix, \code{\link{deSolve}{ode}} method is assumed. The default method overall is \code{"lsodes"} for \code{\link{DAISIE_ML_CS}()} and \code{"ode45"} from \code{\link[deSolve]{ode}()} for @@ -192,7 +194,7 @@ 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 +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. diff --git a/man/DAISIE_ML1.Rd b/man/DAISIE_ML1.Rd index 6386793e..a908afc6 100644 --- a/man/DAISIE_ML1.Rd +++ b/man/DAISIE_ML1.Rd @@ -133,7 +133,7 @@ solvers (steppers) are: \code{"odeint::runge_kutta_fehlberg78"} \code{"odeint::runge_kutta_dopri5"} \code{"odeint::bulirsch_stoer"} -without \code{odeint::}-prefix, \code{\link[deSolve]{ode}} method is +without \code{odeint::}-prefix, \code{\link{deSolve}{ode}} method is assumed. The default method overall is \code{"lsodes"} for \code{\link{DAISIE_ML_CS}()} and \code{"ode45"} from \code{\link[deSolve]{ode}()} for @@ -171,11 +171,13 @@ describing area through time. String checked by \code{\link{is_island_ontogeny_input}()}. \cr In all other functions a numeric describing the type of island ontogeny. Can be \code{0} for constant, \code{1} for a beta function describing area through time. In ML -functions \code{island_ontogeny = NA} assumes constant ontogeny.} +functions \code{island_ontogeny = NA} assumes constant ontogeny. Time +dependent estimation is not yet available as development is still ongoing. +Will return an error if called in that case.} \item{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 +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. diff --git a/man/DAISIE_ML2.Rd b/man/DAISIE_ML2.Rd index e70bb4f2..d7c726c7 100644 --- a/man/DAISIE_ML2.Rd +++ b/man/DAISIE_ML2.Rd @@ -109,7 +109,9 @@ describing area through time. String checked by \code{\link{is_island_ontogeny_input}()}. \cr In all other functions a numeric describing the type of island ontogeny. Can be \code{0} for constant, \code{1} for a beta function describing area through time. In ML -functions \code{island_ontogeny = NA} assumes constant ontogeny.} +functions \code{island_ontogeny = NA} assumes constant ontogeny. Time +dependent estimation is not yet available as development is still ongoing. +Will return an error if called in that case.} \item{tol}{Sets the tolerances in the optimization. Consists of: \cr reltolx = relative tolerance of parameter values in optimization \cr reltolf = @@ -124,7 +126,7 @@ solvers (steppers) are: \code{"odeint::runge_kutta_fehlberg78"} \code{"odeint::runge_kutta_dopri5"} \code{"odeint::bulirsch_stoer"} -without \code{odeint::}-prefix, \code{\link[deSolve]{ode}} method is +without \code{odeint::}-prefix, \code{\link{deSolve}{ode}} method is assumed. The default method overall is \code{"lsodes"} for \code{\link{DAISIE_ML_CS}()} and \code{"ode45"} from \code{\link[deSolve]{ode}()} for @@ -146,7 +148,7 @@ 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 +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. diff --git a/man/DAISIE_ML3.Rd b/man/DAISIE_ML3.Rd deleted file mode 100644 index 20790421..00000000 --- a/man/DAISIE_ML3.Rd +++ /dev/null @@ -1,175 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/DAISIE_ML3.R -\name{DAISIE_ML3} -\alias{DAISIE_ML3} -\title{Computes MLE for single type species under a clade specific scenario with -ontogeny} -\usage{ -DAISIE_ML3( - datalist, - initparsopt, - idparsopt, - parsfix, - idparsfix, - res = 100, - ddmodel = 0, - cond = 0, - island_ontogeny, - tol = c(1e-04, 1e-05, 1e-07), - maxiter = 1000 * round((1.25)^length(idparsopt)), - methode = "lsodes", - optimmethod = "subplex", - CS_version = 1, - verbose = 0, - tolint = c(1e-16, 1e-10), - jitter = 0, - num_cycles = 1 -) -} -\arguments{ -\item{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 the object can of course also be entered directly. -It is an R list object with the following elements.\cr The first element -of the list has two or three components: \cr \cr \code{$island_age} - the -island age \cr Then, depending on whether a distinction between types is -made, we have:\cr \code{$not_present} - the number of mainland lineages -that are not present on the island \cr or:\cr \code{$not_present_type1} - -the number of mainland lineages of type 1 that are not present on the -island \cr \code{$not_present_type2} - the number of mainland lineages of -type 2 that are not present on the island \cr \cr The remaining elements of -the list each contains information on a single colonist lineage on the -island and has 5 components:\cr \cr \code{$colonist_name} - the name of the -species or clade that colonized the island \cr \code{$branching_times} - -island age followed by stem age of the population/species in the case of -Non-endemic, Non-endemic_MaxAge species and Endemic species with no close relatives -on the island. For endemic clades with more than one species on the island -(cladogenetic clades/ radiations) these should be island age followed by the -branching times of the island clade -including the stem age of the clade\cr \code{$stac} - the -status of the colonist \cr \cr * Non_endemic_MaxAge: 1 \cr * Endemic: 2 -\cr * Endemic&Non_Endemic: 3 \cr * Non_Endemic: 4 \cr -* Endemic_Singleton_MaxAge: 5 \cr * Endemic_Clade_MaxAge: 6 -\cr * Endemic&Non_Endemic_Clade_MaxAge: 7 \cr -\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} - -\item{initparsopt}{The initial values of the parameters that must be -optimized, they are all positive.} - -\item{idparsopt}{The ids of the parameters that must be optimized. The ids -are defined as follows: \cr \cr id = 1 corresponds to lambda^c -(cladogenesis rate) \cr id = 2 corresponds to mu (extinction rate) \cr -id = 3 corresponds to K (clade-level carrying capacity) \cr id = 4 -corresponds to gamma (immigration rate) \cr id = 5 corresponds to lambda^a -(anagenesis rate) \cr id = 6 corresponds to lambda^c (cladogenesis rate) -for an optional subset of the species \cr id = 7 corresponds to mu -(extinction rate) for an optional subset of the species\cr id = 8 -corresponds to K (clade-level carrying capacity) for an optional subset of -the species\cr id = 9 corresponds to gamma (immigration rate) for an -optional subset of the species\cr id = 10 corresponds to lambda^a -(anagenesis rate) for an optional subset of the species\cr id = 11 -corresponds to p_f (fraction of mainland species that belongs to the second -subset of species.} - -\item{parsfix}{The values of the parameters that should not be optimized.} - -\item{idparsfix}{The ids of the parameters that should not be optimized, -e.g. c(1,3) if lambda^c and K should not be optimized.} - -\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.} - -\item{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} - -\item{cond}{cond = 0 : conditioning on island age \cr cond = 1 : -conditioning on island age and non-extinction of the island biota \cr. -cond > 1 : conditioning on island age and having at least cond colonizations -on the island. This last option is not yet available for the IW model \cr} - -\item{island_ontogeny}{In \code{\link{DAISIE_sim_time_dep}()}, -\code{\link{DAISIE_ML_CS}} and plotting a string describing the type of -island ontogeny. Can be \code{"const"}, \code{"beta"} for a beta function -describing area through time. String checked by -\code{\link{is_island_ontogeny_input}()}. \cr In all other functions a -numeric describing the type of island ontogeny. Can be \code{0} for -constant, \code{1} for a beta function describing area through time. In ML -functions \code{island_ontogeny = NA} assumes constant ontogeny.} - -\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.} - -\item{maxiter}{Sets the maximum number of iterations in the optimization.} - -\item{methode}{Method of the ODE-solver. Supported Boost \code{ODEINT} -solvers (steppers) are: -\code{"odeint::runge_kutta_cash_karp54"} -\code{"odeint::runge_kutta_fehlberg78"} -\code{"odeint::runge_kutta_dopri5"} -\code{"odeint::bulirsch_stoer"} -without \code{odeint::}-prefix, \code{\link[deSolve]{ode}} method is -assumed. The default method overall is -\code{"lsodes"} for \code{\link{DAISIE_ML_CS}()} -and \code{"ode45"} from \code{\link[deSolve]{ode}()} for -\code{\link{DAISIE_ML_IW}()}.} - -\item{optimmethod}{Method used in likelihood optimization. Default is -`subplex` (see `\link[subplex]{subplex}()` for full details). -Alternative is \code{"simplex"} which was the method in previous versions.} - -\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}{In simulation and dataprep functions a logical, -\code{Default = TRUE} gives intermediate output should be printed. -For ML functions a numeric determining if intermediate output should be -printed, \code{Default = 0} does not print, \code{verbose = 1} prints -intermediate output of the parameters and loglikelihood, \code{verbose = 2} -means also intermediate progress during loglikelihood computation is shown.} - -\item{tolint}{Vector of two elements containing the absolute and relative -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 -maximum loglikelihood. \item{lambda_c}{ gives the maximum likelihood -estimate of lambda^c, the rate of cladogenesis} \item{mu}{ gives the maximum -likelihood estimate of mu, the extinction rate} \item{K}{ gives the maximum -likelihood estimate of K, the carrying-capacity} \item{gamma}{ gives the -maximum likelihood estimate of gamma, the immigration rate } -\item{lambda_a}{ gives the maximum likelihood estimate of lambda^a, the rate -of anagenesis} \item{loglik}{ gives the maximum loglikelihood} -\item{df}{ gives the number of estimated parameters, i.e. degrees of feedom} -\item{conv}{ gives a message on convergence of optimization; -conv = 0 means convergence} -} -\description{ -Computes MLE for single type species under a clade specific scenario with -ontogeny -} -\keyword{internal} diff --git a/man/DAISIE_ML4.Rd b/man/DAISIE_ML4.Rd index fc98fb4f..78c068d8 100644 --- a/man/DAISIE_ML4.Rd +++ b/man/DAISIE_ML4.Rd @@ -107,7 +107,7 @@ solvers (steppers) are: \code{"odeint::runge_kutta_fehlberg78"} \code{"odeint::runge_kutta_dopri5"} \code{"odeint::bulirsch_stoer"} -without \code{odeint::}-prefix, \code{\link[deSolve]{ode}} method is +without \code{odeint::}-prefix, \code{\link{deSolve}{ode}} method is assumed. The default method overall is \code{"lsodes"} for \code{\link{DAISIE_ML_CS}()} and \code{"ode45"} from \code{\link[deSolve]{ode}()} for @@ -145,11 +145,13 @@ describing area through time. String checked by \code{\link{is_island_ontogeny_input}()}. \cr In all other functions a numeric describing the type of island ontogeny. Can be \code{0} for constant, \code{1} for a beta function describing area through time. In ML -functions \code{island_ontogeny = NA} assumes constant ontogeny.} +functions \code{island_ontogeny = NA} assumes constant ontogeny. Time +dependent estimation is not yet available as development is still ongoing. +Will return an error if called in that case.} \item{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 +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. diff --git a/man/DAISIE_ML_IW.Rd b/man/DAISIE_ML_IW.Rd index ff0c6b69..50ef9088 100644 --- a/man/DAISIE_ML_IW.Rd +++ b/man/DAISIE_ML_IW.Rd @@ -105,7 +105,7 @@ solvers (steppers) are: \code{"odeint::runge_kutta_fehlberg78"} \code{"odeint::runge_kutta_dopri5"} \code{"odeint::bulirsch_stoer"} -without \code{odeint::}-prefix, \code{\link[deSolve]{ode}} method is +without \code{odeint::}-prefix, \code{\link{deSolve}{ode}} method is assumed. The default method overall is \code{"lsodes"} for \code{\link{DAISIE_ML_CS}()} and \code{"ode45"} from \code{\link[deSolve]{ode}()} for @@ -127,7 +127,7 @@ 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 +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. diff --git a/man/DAISIE_abm_factor.Rd b/man/DAISIE_abm_factor.Rd new file mode 100644 index 00000000..d11355ed --- /dev/null +++ b/man/DAISIE_abm_factor.Rd @@ -0,0 +1,17 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/DAISIE_loglik_CS.R +\name{DAISIE_abm_factor} +\alias{DAISIE_abm_factor} +\title{Sets or retrieves the factor to calculate the step-size used by the odeint::adams_bashforth[_moulton] solvers.} +\usage{ +DAISIE_abm_factor(factor) +} +\arguments{ +\item{factor}{sets step-size to \code{factor * (t1 - t0)}. \cr} +} +\value{ +current factor +} +\description{ +Sets or retrieves the factor to calculate the step-size used by the odeint::adams_bashforth[_moulton] solvers. +} diff --git a/man/DAISIE_calc_sumstats_pcrates.Rd b/man/DAISIE_calc_sumstats_pcrates.Rd index 69d5a37b..4f8d70d5 100644 --- a/man/DAISIE_calc_sumstats_pcrates.Rd +++ b/man/DAISIE_calc_sumstats_pcrates.Rd @@ -73,7 +73,9 @@ describing area through time. String checked by \code{\link{is_island_ontogeny_input}()}. \cr In all other functions a numeric describing the type of island ontogeny. Can be \code{0} for constant, \code{1} for a beta function describing area through time. In ML -functions \code{island_ontogeny = NA} assumes constant ontogeny.} +functions \code{island_ontogeny = NA} assumes constant ontogeny. Time +dependent estimation is not yet available as development is still ongoing. +Will return an error if called in that case.} \item{sea_level}{In \code{\link{DAISIE_sim_time_dep}()} and plotting a string describing the type of sea level. Can be \code{"const"} or diff --git a/man/DAISIE_loglik_IW.Rd b/man/DAISIE_loglik_IW.Rd index 075c4ab3..45236832 100644 --- a/man/DAISIE_loglik_IW.Rd +++ b/man/DAISIE_loglik_IW.Rd @@ -79,6 +79,8 @@ solvers (steppers) are: \code{'odeint::runge_kutta_fehlberg78'} [default] \code{'odeint::runge_kutta_dopri5'} \code{'odeint::bulirsch_stoer'} +\code{'odeint::adams_bashforth_[1|2|3|4|5|6|7|8]} +\code{'odeint::adams_bashforth_moulton_[1|2|3|4|5|6|7|8]} without \code{odeint::}-prefix, \code{\link[deSolve]{ode}} method is assumed.} diff --git a/man/DAISIE_loglik_integrand.Rd b/man/DAISIE_loglik_integrand.Rd index 98000f50..aec1eaee 100644 --- a/man/DAISIE_loglik_integrand.Rd +++ b/man/DAISIE_loglik_integrand.Rd @@ -51,7 +51,7 @@ solvers (steppers) are: \code{"odeint::runge_kutta_fehlberg78"} \code{"odeint::runge_kutta_dopri5"} \code{"odeint::bulirsch_stoer"} -without \code{odeint::}-prefix, \code{\link[deSolve]{ode}} method is +without \code{odeint::}-prefix, \code{\link{deSolve}{ode}} method is assumed. The default method overall is \code{"lsodes"} for \code{\link{DAISIE_ML_CS}()} and \code{"ode45"} from \code{\link[deSolve]{ode}()} for diff --git a/man/DAISIE_loglik_integrate.Rd b/man/DAISIE_loglik_integrate.Rd index 5b4e64fc..80ba5b8f 100644 --- a/man/DAISIE_loglik_integrate.Rd +++ b/man/DAISIE_loglik_integrate.Rd @@ -56,7 +56,7 @@ solvers (steppers) are: \code{"odeint::runge_kutta_fehlberg78"} \code{"odeint::runge_kutta_dopri5"} \code{"odeint::bulirsch_stoer"} -without \code{odeint::}-prefix, \code{\link[deSolve]{ode}} method is +without \code{odeint::}-prefix, \code{\link{deSolve}{ode}} method is assumed. The default method overall is \code{"lsodes"} for \code{\link{DAISIE_ML_CS}()} and \code{"ode45"} from \code{\link[deSolve]{ode}()} for diff --git a/man/DAISIE_plot_area.Rd b/man/DAISIE_plot_area.Rd index 414672bb..f0330d36 100644 --- a/man/DAISIE_plot_area.Rd +++ b/man/DAISIE_plot_area.Rd @@ -42,7 +42,9 @@ describing area through time. String checked by \code{\link{is_island_ontogeny_input}()}. \cr In all other functions a numeric describing the type of island ontogeny. Can be \code{0} for constant, \code{1} for a beta function describing area through time. In ML -functions \code{island_ontogeny = NA} assumes constant ontogeny.} +functions \code{island_ontogeny = NA} assumes constant ontogeny. Time +dependent estimation is not yet available as development is still ongoing. +Will return an error if called in that case.} \item{resolution}{numeric indicating resolution of plot. Should be < 0.} diff --git a/man/DAISIE_plot_cladogenesis.Rd b/man/DAISIE_plot_cladogenesis.Rd index 5a4fb109..79f95847 100644 --- a/man/DAISIE_plot_cladogenesis.Rd +++ b/man/DAISIE_plot_cladogenesis.Rd @@ -50,7 +50,9 @@ describing area through time. String checked by \code{\link{is_island_ontogeny_input}()}. \cr In all other functions a numeric describing the type of island ontogeny. Can be \code{0} for constant, \code{1} for a beta function describing area through time. In ML -functions \code{island_ontogeny = NA} assumes constant ontogeny.} +functions \code{island_ontogeny = NA} assumes constant ontogeny. Time +dependent estimation is not yet available as development is still ongoing. +Will return an error if called in that case.} \item{sea_level}{In \code{\link{DAISIE_sim_time_dep}()} and plotting a string describing the type of sea level. Can be \code{"const"} or diff --git a/man/DAISIE_plot_extinction.Rd b/man/DAISIE_plot_extinction.Rd index c4b42045..9562e0dd 100644 --- a/man/DAISIE_plot_extinction.Rd +++ b/man/DAISIE_plot_extinction.Rd @@ -56,7 +56,9 @@ describing area through time. String checked by \code{\link{is_island_ontogeny_input}()}. \cr In all other functions a numeric describing the type of island ontogeny. Can be \code{0} for constant, \code{1} for a beta function describing area through time. In ML -functions \code{island_ontogeny = NA} assumes constant ontogeny.} +functions \code{island_ontogeny = NA} assumes constant ontogeny. Time +dependent estimation is not yet available as development is still ongoing. +Will return an error if called in that case.} \item{sea_level}{In \code{\link{DAISIE_sim_time_dep}()} and plotting a string describing the type of sea level. Can be \code{"const"} or diff --git a/man/DAISIE_plot_immigration.Rd b/man/DAISIE_plot_immigration.Rd index 6d8a8ce8..f74e4a57 100644 --- a/man/DAISIE_plot_immigration.Rd +++ b/man/DAISIE_plot_immigration.Rd @@ -65,7 +65,9 @@ describing area through time. String checked by \code{\link{is_island_ontogeny_input}()}. \cr In all other functions a numeric describing the type of island ontogeny. Can be \code{0} for constant, \code{1} for a beta function describing area through time. In ML -functions \code{island_ontogeny = NA} assumes constant ontogeny.} +functions \code{island_ontogeny = NA} assumes constant ontogeny. Time +dependent estimation is not yet available as development is still ongoing. +Will return an error if called in that case.} \item{sea_level}{In \code{\link{DAISIE_sim_time_dep}()} and plotting a string describing the type of sea level. Can be \code{"const"} or diff --git a/man/DAISIE_plot_pc_rates.Rd b/man/DAISIE_plot_pc_rates.Rd index 0cccad65..d7758da5 100644 --- a/man/DAISIE_plot_pc_rates.Rd +++ b/man/DAISIE_plot_pc_rates.Rd @@ -53,7 +53,9 @@ describing area through time. String checked by \code{\link{is_island_ontogeny_input}()}. \cr In all other functions a numeric describing the type of island ontogeny. Can be \code{0} for constant, \code{1} for a beta function describing area through time. In ML -functions \code{island_ontogeny = NA} assumes constant ontogeny.} +functions \code{island_ontogeny = NA} assumes constant ontogeny. Time +dependent estimation is not yet available as development is still ongoing. +Will return an error if called in that case.} \item{sea_level}{In \code{\link{DAISIE_sim_time_dep}()} and plotting a string describing the type of sea level. Can be \code{"const"} or diff --git a/man/DAISIE_sim_core_time_dep.Rd b/man/DAISIE_sim_core_time_dep.Rd index c75c8d0d..f2ba6171 100644 --- a/man/DAISIE_sim_core_time_dep.Rd +++ b/man/DAISIE_sim_core_time_dep.Rd @@ -70,7 +70,9 @@ describing area through time. String checked by \code{\link{is_island_ontogeny_input}()}. \cr In all other functions a numeric describing the type of island ontogeny. Can be \code{0} for constant, \code{1} for a beta function describing area through time. In ML -functions \code{island_ontogeny = NA} assumes constant ontogeny.} +functions \code{island_ontogeny = NA} assumes constant ontogeny. Time +dependent estimation is not yet available as development is still ongoing. +Will return an error if called in that case.} \item{sea_level}{In \code{\link{DAISIE_sim_time_dep}()} and plotting a string describing the type of sea level. Can be \code{"const"} or diff --git a/man/DAISIE_sim_core_trait_dep.Rd b/man/DAISIE_sim_core_trait_dep.Rd index bdaa60e8..36faea6f 100644 --- a/man/DAISIE_sim_core_trait_dep.Rd +++ b/man/DAISIE_sim_core_trait_dep.Rd @@ -60,7 +60,9 @@ describing area through time. String checked by \code{\link{is_island_ontogeny_input}()}. \cr In all other functions a numeric describing the type of island ontogeny. Can be \code{0} for constant, \code{1} for a beta function describing area through time. In ML -functions \code{island_ontogeny = NA} assumes constant ontogeny.} +functions \code{island_ontogeny = NA} assumes constant ontogeny. Time +dependent estimation is not yet available as development is still ongoing. +Will return an error if called in that case.} \item{sea_level}{In \code{\link{DAISIE_sim_time_dep}()} and plotting a string describing the type of sea level. Can be \code{"const"} or diff --git a/man/DAISIE_sim_time_dep.Rd b/man/DAISIE_sim_time_dep.Rd index e8e61662..00f29ee9 100644 --- a/man/DAISIE_sim_time_dep.Rd +++ b/man/DAISIE_sim_time_dep.Rd @@ -119,7 +119,9 @@ describing area through time. String checked by \code{\link{is_island_ontogeny_input}()}. \cr In all other functions a numeric describing the type of island ontogeny. Can be \code{0} for constant, \code{1} for a beta function describing area through time. In ML -functions \code{island_ontogeny = NA} assumes constant ontogeny.} +functions \code{island_ontogeny = NA} assumes constant ontogeny. Time +dependent estimation is not yet available as development is still ongoing. +Will return an error if called in that case.} \item{sea_level}{In \code{\link{DAISIE_sim_time_dep}()} and plotting a string describing the type of sea level. Can be \code{"const"} or diff --git a/man/DAISIE_sim_time_dep_cs.Rd b/man/DAISIE_sim_time_dep_cs.Rd index 32fc911f..66b8180c 100644 --- a/man/DAISIE_sim_time_dep_cs.Rd +++ b/man/DAISIE_sim_time_dep_cs.Rd @@ -98,7 +98,9 @@ describing area through time. String checked by \code{\link{is_island_ontogeny_input}()}. \cr In all other functions a numeric describing the type of island ontogeny. Can be \code{0} for constant, \code{1} for a beta function describing area through time. In ML -functions \code{island_ontogeny = NA} assumes constant ontogeny.} +functions \code{island_ontogeny = NA} assumes constant ontogeny. Time +dependent estimation is not yet available as development is still ongoing. +Will return an error if called in that case.} \item{sea_level}{In \code{\link{DAISIE_sim_time_dep}()} and plotting a string describing the type of sea level. Can be \code{"const"} or diff --git a/man/DAISIE_sim_time_dep_gw.Rd b/man/DAISIE_sim_time_dep_gw.Rd index 7c3e5e64..9afbf070 100644 --- a/man/DAISIE_sim_time_dep_gw.Rd +++ b/man/DAISIE_sim_time_dep_gw.Rd @@ -102,7 +102,9 @@ describing area through time. String checked by \code{\link{is_island_ontogeny_input}()}. \cr In all other functions a numeric describing the type of island ontogeny. Can be \code{0} for constant, \code{1} for a beta function describing area through time. In ML -functions \code{island_ontogeny = NA} assumes constant ontogeny.} +functions \code{island_ontogeny = NA} assumes constant ontogeny. Time +dependent estimation is not yet available as development is still ongoing. +Will return an error if called in that case.} \item{sea_level}{In \code{\link{DAISIE_sim_time_dep}()} and plotting a string describing the type of sea level. Can be \code{"const"} or diff --git a/man/DAISIE_sim_time_dep_iw.Rd b/man/DAISIE_sim_time_dep_iw.Rd index 90c5b4d0..602576e8 100644 --- a/man/DAISIE_sim_time_dep_iw.Rd +++ b/man/DAISIE_sim_time_dep_iw.Rd @@ -98,7 +98,9 @@ describing area through time. String checked by \code{\link{is_island_ontogeny_input}()}. \cr In all other functions a numeric describing the type of island ontogeny. Can be \code{0} for constant, \code{1} for a beta function describing area through time. In ML -functions \code{island_ontogeny = NA} assumes constant ontogeny.} +functions \code{island_ontogeny = NA} assumes constant ontogeny. Time +dependent estimation is not yet available as development is still ongoing. +Will return an error if called in that case.} \item{sea_level}{In \code{\link{DAISIE_sim_time_dep}()} and plotting a string describing the type of sea level. Can be \code{"const"} or diff --git a/man/DAISIE_sim_trait_dep.Rd b/man/DAISIE_sim_trait_dep.Rd index 172bdfb7..fada74dc 100644 --- a/man/DAISIE_sim_trait_dep.Rd +++ b/man/DAISIE_sim_trait_dep.Rd @@ -92,7 +92,9 @@ describing area through time. String checked by \code{\link{is_island_ontogeny_input}()}. \cr In all other functions a numeric describing the type of island ontogeny. Can be \code{0} for constant, \code{1} for a beta function describing area through time. In ML -functions \code{island_ontogeny = NA} assumes constant ontogeny.} +functions \code{island_ontogeny = NA} assumes constant ontogeny. Time +dependent estimation is not yet available as development is still ongoing. +Will return an error if called in that case.} \item{sea_level}{In \code{\link{DAISIE_sim_time_dep}()} and plotting a string describing the type of sea level. Can be \code{"const"} or diff --git a/man/add_column_to_dataframe.Rd b/man/add_column_to_dataframe.Rd new file mode 100644 index 00000000..5e6cfc35 --- /dev/null +++ b/man/add_column_to_dataframe.Rd @@ -0,0 +1,24 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/DAISIE_utils.R +\name{add_column_to_dataframe} +\alias{add_column_to_dataframe} +\title{Add a column to a data frame} +\usage{ +add_column_to_dataframe(df, position, column_to_insert) +} +\arguments{ +\item{df}{data frame to add the column to} + +\item{position}{location in data frame where to insert the column. +Position can also be a name of a column} + +\item{column_to_insert}{the elements of the column to insert. If +the column has a name, this name will be copied into the data frame. +Id is does not have a name, it will get the name "nc".} +} +\value{ +A data frame with the column inserted +} +\description{ +Add a column to a data frame +} diff --git a/man/default_params_doc.Rd b/man/default_params_doc.Rd index ceda5ae4..1cd4d25c 100644 --- a/man/default_params_doc.Rd +++ b/man/default_params_doc.Rd @@ -254,7 +254,9 @@ describing area through time. String checked by \code{\link{is_island_ontogeny_input}()}. \cr In all other functions a numeric describing the type of island ontogeny. Can be \code{0} for constant, \code{1} for a beta function describing area through time. In ML -functions \code{island_ontogeny = NA} assumes constant ontogeny.} +functions \code{island_ontogeny = NA} assumes constant ontogeny. Time +dependent estimation is not yet available as development is still ongoing. +Will return an error if called in that case.} \item{sea_level}{In \code{\link{DAISIE_sim_time_dep}()} and plotting a string describing the type of sea level. Can be \code{"const"} or @@ -461,7 +463,7 @@ solvers (steppers) are: \code{"odeint::runge_kutta_fehlberg78"} \code{"odeint::runge_kutta_dopri5"} \code{"odeint::bulirsch_stoer"} -without \code{odeint::}-prefix, \code{\link[deSolve]{ode}} method is +without \code{odeint::}-prefix, \code{\link{deSolve}{ode}} method is assumed. The default method overall is \code{"lsodes"} for \code{\link{DAISIE_ML_CS}()} and \code{"ode45"} from \code{\link[deSolve]{ode}()} for @@ -665,7 +667,7 @@ at the end of the simulation).} \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 +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. diff --git a/man/get_clado_rate_per_capita.Rd b/man/get_clado_rate_per_capita.Rd new file mode 100644 index 00000000..66c08deb --- /dev/null +++ b/man/get_clado_rate_per_capita.Rd @@ -0,0 +1,38 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/DAISIE_rates.R +\name{get_clado_rate_per_capita} +\alias{get_clado_rate_per_capita} +\title{Calculate per-capita cladogenesis rate} +\usage{ +get_clado_rate_per_capita(lac, d, num_spec, K, A = 1) +} +\arguments{ +\item{lac}{A numeric with the per capita cladogenesis rate.} + +\item{d}{Numeric defining the scaling parameter for exponent for +calculating cladogenesis rate.} + +\item{num_spec}{A numeric with the current number of species.} + +\item{K}{A numeric with carrying capacity.} + +\item{A}{A numeric value for island area at a given point in time.} +} +\value{ +Numeric with the per-capita cladogenesis rate given a base +cladogenesis rate, K, A and the d hyperparameter. +} +\description{ +This function is only called directly inside the RHS of the ontogeny +likelihood functions. In all other cases \code{\link{get_clado_rate}()} is to +be called instead. +} +\examples{ +lac <- 0.4 +d <- 0 +num_spec <- 2 +K <- 10 +A <- 1 +clado_rate_pc <- DAISIE:::get_clado_rate_per_capita(lac, d, num_spec, K, A) +} +\keyword{internal} diff --git a/man/get_ext_rate.Rd b/man/get_ext_rate.Rd index 79a5c6de..ca6e2dbc 100644 --- a/man/get_ext_rate.Rd +++ b/man/get_ext_rate.Rd @@ -2,14 +2,14 @@ % Please edit documentation in R/DAISIE_rates.R \name{get_ext_rate} \alias{get_ext_rate} -\title{Function to describe changes in extinction rate through time.} +\title{Calculate extinction rate} \usage{ get_ext_rate( mu, hyper_pars, extcutoff = 1000, num_spec, - A, + A = 1, trait_pars = NULL, island_spec = NULL ) @@ -48,8 +48,13 @@ rate preventing it from being too large and slowing down simulation.} \item{island_spec}{Matrix with current state of simulation containing number of species.} } +\value{ +A numeric, with the extinction rate given the base extinction rate, +if present, the hyperparemeter \code{x}, A(t) if time-dependent and traits +if running a traint model. +} \description{ -Function to describe changes in extinction rate through time. +Calculate extinction rate } \references{ Valente, Luis M., Rampal S. Etienne, and Albert B. Phillimore. diff --git a/man/get_ext_rate_per_capita.Rd b/man/get_ext_rate_per_capita.Rd new file mode 100644 index 00000000..66018299 --- /dev/null +++ b/man/get_ext_rate_per_capita.Rd @@ -0,0 +1,34 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/DAISIE_rates.R +\name{get_ext_rate_per_capita} +\alias{get_ext_rate_per_capita} +\title{Function to describe per-capita changes in extinction rate through time} +\usage{ +get_ext_rate_per_capita(mu, x, extcutoff = 1000, A = 1) +} +\arguments{ +\item{mu}{A numeric with the per capita extinction rate.} + +\item{x}{Numeric defining the exponent for calculating extinction rate.} + +\item{extcutoff}{A numeric with the cutoff for the the maximum extinction +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.} +} +\value{ +Numeric with per capita extinction rate, given A(t), x, and mu0. +} +\description{ +This function is only called directly inside the RHS of the ontogeny +likelihood functions. In all other cases \code{\link{get_ext_rate}()} is to +be called instead. +} +\examples{ +ext_rate_per_capita <- DAISIE:::get_ext_rate_per_capita( + mu = 0.5, + x = 1, + A = 1000 +) +} +\keyword{internal} diff --git a/man/get_global_max_area.Rd b/man/get_global_max_area.Rd index ecc55c55..76685f40 100644 --- a/man/get_global_max_area.Rd +++ b/man/get_global_max_area.Rd @@ -35,7 +35,9 @@ describing area through time. String checked by \code{\link{is_island_ontogeny_input}()}. \cr In all other functions a numeric describing the type of island ontogeny. Can be \code{0} for constant, \code{1} for a beta function describing area through time. In ML -functions \code{island_ontogeny = NA} assumes constant ontogeny.} +functions \code{island_ontogeny = NA} assumes constant ontogeny. Time +dependent estimation is not yet available as development is still ongoing. +Will return an error if called in that case.} \item{sea_level}{In \code{\link{DAISIE_sim_time_dep}()} and plotting a string describing the type of sea level. Can be \code{"const"} or diff --git a/man/get_global_min_area.Rd b/man/get_global_min_area.Rd index 1e63d972..bfb80a7d 100644 --- a/man/get_global_min_area.Rd +++ b/man/get_global_min_area.Rd @@ -35,7 +35,9 @@ describing area through time. String checked by \code{\link{is_island_ontogeny_input}()}. \cr In all other functions a numeric describing the type of island ontogeny. Can be \code{0} for constant, \code{1} for a beta function describing area through time. In ML -functions \code{island_ontogeny = NA} assumes constant ontogeny.} +functions \code{island_ontogeny = NA} assumes constant ontogeny. Time +dependent estimation is not yet available as development is still ongoing. +Will return an error if called in that case.} \item{sea_level}{In \code{\link{DAISIE_sim_time_dep}()} and plotting a string describing the type of sea level. Can be \code{"const"} or diff --git a/man/get_immig_rate.Rd b/man/get_immig_rate.Rd index d367f752..241b6dca 100644 --- a/man/get_immig_rate.Rd +++ b/man/get_immig_rate.Rd @@ -6,7 +6,7 @@ \usage{ get_immig_rate( gam, - A, + A = 1, num_spec, K, mainland_n, diff --git a/man/get_immig_rate_per_capita.Rd b/man/get_immig_rate_per_capita.Rd new file mode 100644 index 00000000..2340f481 --- /dev/null +++ b/man/get_immig_rate_per_capita.Rd @@ -0,0 +1,34 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/DAISIE_rates.R +\name{get_immig_rate_per_capita} +\alias{get_immig_rate_per_capita} +\title{Calculate per-capita immigration rate} +\usage{ +get_immig_rate_per_capita(gam, num_spec, K, A = 1) +} +\arguments{ +\item{gam}{A numeric with the per capita immigration rate.} + +\item{num_spec}{A numeric with the current number of species.} + +\item{K}{A numeric with carrying capacity.} + +\item{A}{A numeric value for island area at a given point in time.} +} +\value{ +A numeric with the per-capita immigration rate given A(t) and K. +} +\description{ +This function is only called directly inside the RHS of the ontogeny +likelihood functions. In all other cases \code{\link{get_immig_rate}()} is to +be called instead. +} +\examples{ +immig_rate_per_capita <- DAISIE:::get_immig_rate_per_capita( + gam = 0.001, + num_spec = 5, + K = 20, + A = 1000 +) +} +\keyword{internal} diff --git a/man/is_island_ontogeny_input.Rd b/man/is_island_ontogeny_input.Rd index 471f5192..847b6877 100644 --- a/man/is_island_ontogeny_input.Rd +++ b/man/is_island_ontogeny_input.Rd @@ -14,7 +14,9 @@ describing area through time. String checked by \code{\link{is_island_ontogeny_input}()}. \cr In all other functions a numeric describing the type of island ontogeny. Can be \code{0} for constant, \code{1} for a beta function describing area through time. In ML -functions \code{island_ontogeny = NA} assumes constant ontogeny.} +functions \code{island_ontogeny = NA} assumes constant ontogeny. Time +dependent estimation is not yet available as development is still ongoing. +Will return an error if called in that case.} } \value{ Boolean stating if island_ontogeny is correct. diff --git a/man/island_area.Rd b/man/island_area.Rd index 79e29c80..d0ca12f9 100644 --- a/man/island_area.Rd +++ b/man/island_area.Rd @@ -37,7 +37,9 @@ describing area through time. String checked by \code{\link{is_island_ontogeny_input}()}. \cr In all other functions a numeric describing the type of island ontogeny. Can be \code{0} for constant, \code{1} for a beta function describing area through time. In ML -functions \code{island_ontogeny = NA} assumes constant ontogeny.} +functions \code{island_ontogeny = NA} assumes constant ontogeny. Time +dependent estimation is not yet available as development is still ongoing. +Will return an error if called in that case.} \item{sea_level}{In \code{\link{DAISIE_sim_time_dep}()} and plotting a string describing the type of sea level. Can be \code{"const"} or diff --git a/man/island_area_vector.Rd b/man/island_area_vector.Rd index d2b46ae6..9ec9a62a 100644 --- a/man/island_area_vector.Rd +++ b/man/island_area_vector.Rd @@ -4,7 +4,14 @@ \alias{island_area_vector} \title{Computes island_area, but takes vector as argument (needed by )} \usage{ -island_area_vector(timeval, area_pars, island_ontogeny, sea_level) +island_area_vector( + timeval, + area_pars, + island_ontogeny, + sea_level, + total_time, + peak +) } \arguments{ \item{timeval}{current time of simulation} diff --git a/man/translate_island_ontogeny.Rd b/man/translate_island_ontogeny.Rd index 9c4091f0..b491a2e6 100644 --- a/man/translate_island_ontogeny.Rd +++ b/man/translate_island_ontogeny.Rd @@ -14,7 +14,9 @@ describing area through time. String checked by \code{\link{is_island_ontogeny_input}()}. \cr In all other functions a numeric describing the type of island ontogeny. Can be \code{0} for constant, \code{1} for a beta function describing area through time. In ML -functions \code{island_ontogeny = NA} assumes constant ontogeny.} +functions \code{island_ontogeny = NA} assumes constant ontogeny. Time +dependent estimation is not yet available as development is still ongoing. +Will return an error if called in that case.} } \value{ Numeric, 0 for null-ontogeny, 1 for beta function diff --git a/man/update_rates.Rd b/man/update_rates.Rd index eac7a1aa..0e9c3b2f 100644 --- a/man/update_rates.Rd +++ b/man/update_rates.Rd @@ -72,7 +72,9 @@ describing area through time. String checked by \code{\link{is_island_ontogeny_input}()}. \cr In all other functions a numeric describing the type of island ontogeny. Can be \code{0} for constant, \code{1} for a beta function describing area through time. In ML -functions \code{island_ontogeny = NA} assumes constant ontogeny.} +functions \code{island_ontogeny = NA} assumes constant ontogeny. Time +dependent estimation is not yet available as development is still ongoing. +Will return an error if called in that case.} \item{sea_level}{In \code{\link{DAISIE_sim_time_dep}()} and plotting a string describing the type of sea level. Can be \code{"const"} or diff --git a/src/DAISIE_CS.cpp b/src/DAISIE_CS.cpp index dfce9ef4..902ffdeb 100644 --- a/src/DAISIE_CS.cpp +++ b/src/DAISIE_CS.cpp @@ -17,6 +17,9 @@ namespace { static constexpr int default_max_cs_steps = 1'000'000; static int max_cs_steps = default_max_cs_steps; + // step-size factor for adams_bashforth_moulton integration + static constexpr double default_abm_factor = 0.0001; + static double abm_factor = default_abm_factor; // class padded_vector_view @@ -73,8 +76,8 @@ namespace { // xx2 = c(0,0,x[(lx + 1):(2 * lx)],0) // xx3 = x[2 * lx + 1] // using padded views instead of vectors: - const auto xx1 = padded_vector_view(2, x.data(), p_.lx); - const auto xx2 = padded_vector_view(2, x.data() + p_.lx, p_.lx); + const auto xx1 = padded_vector_view(2, x.data().begin(), p_.lx); + const auto xx2 = padded_vector_view(2, x.data().begin() + p_.lx, p_.lx); const auto xx3 = x[2 * p_.lx]; // DO I = 1, N + 4 + 2 * kk @@ -86,11 +89,11 @@ namespace { // ENDDO // using views instead of vectors: const auto chunk = p_.lx + 4 + 2 * p_.kk; - const auto laavec = p_.P.data(); - const auto lacvec = p_.P.data() + chunk; - const auto muvec = p_.P.data() + 2 * chunk; - const auto gamvec = p_.P.data() + 3 * chunk; - const auto nn = p_.P.data() + 4 * chunk; + const auto laavec = p_.P.data().begin(); + const auto lacvec = p_.P.data().begin() + chunk; + const auto muvec = p_.P.data().begin() + 2 * chunk; + const auto gamvec = p_.P.data().begin() + 3 * chunk; + const auto nn = p_.P.data().begin() + 4 * chunk; // DO I = 3, N + 2 // il1(I - 2) = I + kk - 1 @@ -133,7 +136,7 @@ namespace { // dConc(N + I) = FF1 // ENDDO // using views into output vector: - auto dx1 = dx.data(); + auto dx1 = dx.data().begin(); auto dx2 = dx1 + p_.lx; for (int i = 0; i < p_.lx; ++i) { dx1[i] = laavec[il1 + i + 1] * xx2[ix1 + i] @@ -154,9 +157,147 @@ namespace { // dConc(1) = dConc(1) + laavec(il3in3(1)) * xx3 // dConc(2) = dConc(2) + 2 * lacvec(il3in3(1)) * xx3 // ENDIF - if (1 == p_.kk) { - dx1[0] += laavec[il3in3] * xx3; - dx2[1] += 2.0 * lacvec[il3in3] * xx3; + // if (1 == p_.kk) { + // dx1[0] += laavec[il3in3] * xx3; + // dx2[1] += 2.0 * lacvec[il3in3] * xx3; + //} + + // FFF = laavec(il3in3(1)) + lacvec(il3in3(1)) + // FFF = FFF + gamvec(il3in3(1)) + muvec(il3in3(1)) + // dConc(2 * N + 1) = -1 * FFF * xx3 + auto dx3 = dx2 + p_.lx; + dx3[0] = -(laavec[il3in3] + lacvec[il3in3] + gamvec[il3in3] + muvec[il3in3]) * xx3; + } + + private: + const param_t p_; + }; + + class cpp_daisie_cs_runmod_1 + { + public: + cpp_daisie_cs_runmod_1(param_t&& p) : + p_(p) + { + } + + // odeint interface + void operator()(const state_type& x, state_type& dx, double) const + { + if (++p_.steps > max_cs_steps) throw std::runtime_error("cpp_daisie_cs_runmod_1: too many steps"); + + // xx1 = c(0,0,x[1:lx],0) + // xx2 = c(0,0,x[(lx + 1):(2 * lx)],0) + // xx3 = c(0,0,x[(2 * lx + 1):(3 * lx)],0) + // xx4 <- c(0,0,x[(3 * lx + 1):(4 * lx)],0) + // using padded views instead of vectors: + const auto xx1 = padded_vector_view(2, x.data().begin(), p_.lx); + const auto xx2 = padded_vector_view(2, x.data().begin() + p_.lx, p_.lx); + const auto xx3 = padded_vector_view(2, x.data().begin() + 2 * p_.lx, p_.lx); + const auto xx4 = padded_vector_view(2, x.data().begin() + 3 * p_.lx, p_.lx); + + // DO I = 1, N + 4 + 2 * kk + // laavec(I) = P(I) + // lacvec(I) = P(I + N + 4 + 2 * kk) + // muvec(I) = P(I + 2 * (N + 4 + 2 * kk)) + // gamvec(I) = P(I + 3 * (N + 4 + 2 * kk)) + // nn(I) = P(I + 4 * (N + 4 + 2 * kk)) + // ENDDO + // using views instead of vectors: + const auto chunk = p_.lx + 4 + 2 * p_.kk; + const auto laavec = p_.P.data().begin(); + const auto lacvec = p_.P.data().begin() + chunk; + const auto muvec = p_.P.data().begin() + 2 * chunk; + const auto gamvec = p_.P.data().begin() + 3 * chunk; + const auto nn = p_.P.data().begin() + 4 * chunk; + + // DO I = 3, N + 2 + // il1(I - 2) = I + kk - 1 + // il2(I - 2) = I + kk + 1 + // il3in3(I - 2) = I + kk + // il4(I - 2) = I + kk - 2 + // in1(I - 2) = I + 2 * kk - 1 + // in2ix2(I - 2) = I + 1 + // in4ix1(I - 2) = I - 1 + // ix3(I - 2) = I + // ix4(I - 2) = I - 2 + // ENDDO + // using offsets into our views instead of vectors: + const int il1 = 2 + p_.kk - 1; + const int il2 = 2 + p_.kk + 1; + const int il3in3 = 2 + p_.kk; + const int il4 = 2 + p_.kk - 2; + const int in1 = 2 + 2 * p_.kk - 1; + const int in2ix2 = 2 + 1; // spilt in in2, ix2 at no cost? + const int in4ix1 = 2 - 1; // split in in4, ix1 at no cost? + const int ix3 = 2; + const int ix4 = 2 - 2; + + //DO I = 1, N + // FF1 = lacvec(il1(I)) * xx1(in4ix1(I)) + // FF1 = FF1 + laavec(il1(I) + 1) * xx2(in4ix1(I)) + // FF1 = FF1 + lacvec(il4(I) + 1) * xx2(ix4(I)) + // FF1 = FF1 + muvec(il2(I)) * nn(in2ix2(I)) * xx1(in2ix2(I)) + // FF1 = FF1 + muvec(il3in3(I) + 1) * xx2(ix3(I)) + // FFF = muvec(il3in3(I)) + lacvec(il3in3(I)) + // FF1 = FF1 - FFF * nn(il3in3(I)) * xx1(ix3(I)) + // dConc(I) = FF1 - gamvec(il3in3(I)) * xx1(ix3(I)) + // FF1 = gamvec(il3in3(I)) * xx1(ix3(I)) + // FF1 = FF1 + gamvec(il3in3(I)) * xx3(ix3(I)) + // FF1 = FF1 + gamvec(il3in3(I) + 1) * xx4(ix3(I)) + // FF1 = FF1 + lacvec(il1(I) + 1) * nn(in1(I)) * xx2(in4ix1(I)) + // FF1 = FF1 + muvec(il2(I) + 1) * nn(in2ix2(I)) * xx2(in2ix2(I)) + // FFF = muvec(il3in3(I) + 1) + lacvec(il3in3(I) + 1) + // FF1 = FF1 - FFF * nn(il3in3(I) + 1) * xx2(ix3(I)) + // FF1 = FF1 - laavec(il3in3(I) + 1) * xx2(ix3(I)) + // dConc(N + I) = FF1 + // FF1 = lacvec(il1(I)) * nn(in1(I)) * xx3(in4ix1(I)) + // FF1 = FF1 + laavec(il1(I) + 1) * xx4(in4ix1(I)) + // FF1 = FF1 + lacvec(il4(I) + 1) * xx4(ix4(I)) + // FF1 = FF1 + muvec(il2(I)) * nn(in2ix2(I)) * xx3(in2ix2(I)) + // FF1 = FF1 + muvec(il3in3(I) + 1) * xx4(ix3(I)) + // FFF = lacvec(il3in3(I)) + muvec(il3in3(I)) + // FF1 = FF1 - FFF * nn(il3in3(I)) * xx3(ix3(I)) + // FF1 = FF1 - gamvec(il3in3(I)) * xx3(ix3(I)) + // dConc(2 * N + I) = FF1 + // FF1 = lacvec(il1(I) + 1) * nn(in1(I)) * xx4(in4ix1(I)) + // FF1 = FF1 + muvec(il2(I) + 1) * nn(in2ix2(I)) * xx4(in2ix2(I)) + // FFF = lacvec(il3in3(I) + 1) + muvec(il3in3(I) + 1) + // FF1 = FF1 - FFF * nn(il3in3(I) + 1) * xx4(ix3(I)) + // FF1 = FF1 - gamvec(il3in3(I) + 1) * xx4(ix3(I)) + // dConc(3 * N + I) = FF1 + //ENDDO + // using views into output vector: + auto dx1 = dx.data().begin(); + auto dx2 = dx1 + p_.lx; + auto dx3 = dx2 + p_.lx; + auto dx4 = dx3 + p_.lx; + for (int i = 0; i < p_.lx; ++i) { + dx1[i] = lacvec[il1 + i] * nn[in1 + i] * xx1[in4ix1 + i] + + laavec[il1 + i + 1] * xx2[in4ix1 + i] + + lacvec[il4 + i + 1] * xx2[ix4 + i] + + muvec[il2 + i] * nn[in2ix2 + i] * xx1[in2ix2 + i] + + muvec[il3in3 + i + 1] * xx2[ix3 + i] + - (muvec[il3in3 + i] + lacvec[il3in3 + i]) * nn[il3in3 + i] * xx1[ix3 + i] + - gamvec[il3in3 + i] * xx1[ix3 + i]; + dx2[i] = gamvec[il3in3 + i] * xx1[ix3 + i] + + gamvec[il3in3 + i] * xx3[ix3 + i] + + gamvec[il3in3 + i + 1] * xx4[ix3 + i] + + lacvec[il1 + i + 1] * nn[in1 + i] * xx2[in4ix1 + i] + + muvec[il2 + i + 1] * nn[in2ix2 + i] * xx2[in2ix2 + i] + - (muvec[il3in3 + 1 + i] + lacvec[il3in3 + 1 + i]) * nn[il3in3 + i + 1] * xx2[ix3 + i] + - laavec[il3in3 + i] * xx2[ix3 + i]; + dx3[i] = lacvec[il1 + i] * nn[in1 + i] * xx3[in4ix1 + i] + + laavec[il1 + i + 1] * xx4[in4ix1 + i] + + lacvec[il4 + i + 1] * xx4[ix4 + i] + + muvec[il2 + i] * nn[in2ix2 + i] * xx3[in2ix2 + i] + + muvec[il3in3 + i + 1] * xx4[ix3 + i] + - (lacvec[il3in3 + i] + muvec[il3in3 + i]) * nn[il3in3 + i] * xx3[ix3 + i] + - gamvec[il3in3 + i] * xx3[ix3 + i]; + dx4[i] = lacvec[il1 + i + 1] * nn[in1 + i] * xx4[in4ix1 + i] + + muvec[il2 + i + 1] * nn[in2ix2 + i] * xx4[in2ix2 + i] + - (lacvec[il3in3 + i + 1] + muvec[il3in3 + i + 1]) * nn[il3in3 + i + 1] * xx4[ix3 + i] + - gamvec[il3in3 + i + 1] * xx4[ix3 + i]; } } @@ -182,9 +323,9 @@ namespace { // xx2 = c(0,0,x[(lx + 1):(2 * lx)],0) // xx3 = c(0,0,x[(2 * lx + 1):(3 * lx)],0) // using padded views instead of vectors: - const auto xx1 = padded_vector_view(2, x.data(), p_.lx); - const auto xx2 = padded_vector_view(2, x.data() + p_.lx, p_.lx); - const auto xx3 = padded_vector_view(2, x.data() + 2 * p_.lx, p_.lx); + const auto xx1 = padded_vector_view(2, x.data().begin(), p_.lx); + const auto xx2 = padded_vector_view(2, x.data().begin() + p_.lx, p_.lx); + const auto xx3 = padded_vector_view(2, x.data().begin() + 2 * p_.lx, p_.lx); // DO I = 1, N + 4 + 2 * kk // laavec(I) = P(I) @@ -195,11 +336,11 @@ namespace { // ENDDO // using views instead of vectors: const auto chunk = p_.lx + 4 + 2 * p_.kk; - const auto laavec = p_.P.data(); - const auto lacvec = p_.P.data() + chunk; - const auto muvec = p_.P.data() + 2 * chunk; - const auto gamvec = p_.P.data() + 3 * chunk; - const auto nn = p_.P.data() + 4 * chunk; + const auto laavec = p_.P.data().begin(); + const auto lacvec = p_.P.data().begin() + chunk; + const auto muvec = p_.P.data().begin() + 2 * chunk; + const auto gamvec = p_.P.data().begin() + 3 * chunk; + const auto nn = p_.P.data().begin() + 4 * chunk; // DO I = 3, N + 2 // il1(I - 2) = I + kk - 1 @@ -253,7 +394,7 @@ namespace { // dConc(2 * N + I) = FF1 // ENDDO // using views into output vector: - auto dx1 = dx.data(); + auto dx1 = dx.data().begin(); auto dx2 = dx1 + p_.lx; auto dx3 = dx2 + p_.lx; for (int i = 0; i < p_.lx; ++i) { @@ -266,7 +407,6 @@ namespace { - gamvec[il3in3 + i] * xx1[ix3 + i]; if (1 == p_.kk) { dx1[i] += laavec[il3in3 + i] * xx3[ix3 + i] + 2.0 * lacvec[il1 + i] * xx3[in4ix1 + i]; - } dx2[i] = gamvec[il3in3 + i] * xx1[ix3 + i] + lacvec[il1 + i + 1] * nn[in1 + i] * xx2[in4ix1 + i] @@ -295,7 +435,7 @@ BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; auto runmod = as(rrunmod); - auto y = as>(ry); + auto y = as(ry); auto times = as>(rtimes); auto lx = as(rlx); auto kk = as(rkk); @@ -308,6 +448,10 @@ BEGIN_RCPP cpp_daisie_cs_runmod rhs(std::move(p)); daisie_odeint::integrate(stepper, std::ref(rhs), y, times[0], times[1], atol, rtol); } + else if (runmod == "daisie_runmod1") { + cpp_daisie_cs_runmod_1 rhs(std::move(p)); + daisie_odeint::integrate(stepper, std::ref(rhs), y, times[0], times[1], atol, rtol); + } else if (runmod == "daisie_runmod2") { cpp_daisie_cs_runmod_2 rhs(std::move(p)); daisie_odeint::integrate(stepper, std::ref(rhs), y, times[0], times[1], atol, rtol); @@ -329,3 +473,21 @@ RcppExport SEXP daisie_odeint_cs_max_steps(SEXP rmax_steps) { END_RCPP } + +namespace daisie_odeint { + + // step-size factor for adams_bashforth_moulton integration + constexpr double default_abm_factor = 0.0001; + double abm_factor = default_abm_factor; + +} + + +// misplaced +RcppExport SEXP daisie_odeint_abm_factor(SEXP rfactor) { + BEGIN_RCPP + daisie_odeint::abm_factor = (0 < as(rfactor)) ? as(rfactor) : daisie_odeint::default_abm_factor; + return wrap(daisie_odeint::abm_factor); + END_RCPP +} + diff --git a/src/DAISIE_IW.cpp b/src/DAISIE_IW.cpp index f6ef071f..e6f55b53 100644 --- a/src/DAISIE_IW.cpp +++ b/src/DAISIE_IW.cpp @@ -1,11 +1,7 @@ -// [[Rcpp::plugins(cpp14)]] -// [[Rcpp::plugins(openmp)]] -// [[Rcpp::depends(RcppEigen)]] -// [[Rcpp::depends(BH)]] - - //' @export daisie_odeint_iw +// [[Rcpp::plugins(openmp)]] +// [[Rcpp::depends(RcppEigen)]] #include "DAISIE_odeint.h" #define EIGEN_USE_THREADS @@ -206,8 +202,8 @@ namespace { // odeint interface void operator()(const state_type& x, state_type& dxdt, double) { - (iw2) ? iw2->rhs(x.data(), dxdt.data(), dev.get()) - : iw3->rhs(x.data(), dxdt.data(), dev.get()); + (iw2) ? iw2->rhs(x.data().begin(), dxdt.data().begin(), dev.get()) + : iw3->rhs(x.data().begin(), dxdt.data().begin(), dev.get()); } }; @@ -222,7 +218,7 @@ BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; auto y = as(ry); - auto times = as(rtimes); + auto times = as>(rtimes); auto pars = as(rpars); auto stepper = as(Stepper); auto atol = as(atolint); diff --git a/src/DAISIE_loglik_rhs_FORTRAN.f95 b/src/DAISIE_loglik_rhs_FORTRAN.f95 index bf449ebe..2bb49994 100644 --- a/src/DAISIE_loglik_rhs_FORTRAN.f95 +++ b/src/DAISIE_loglik_rhs_FORTRAN.f95 @@ -1,12 +1,5 @@ - -! ------------------------------------------------------------------------ -! ------------------------------------------------------------------------ -! Example how to apply Fortran code with variable-length parameters -! --------------------------------------------------------------------- -! ------------------------------------------------------------------------ - !========================================================================== -! Helper function: +! Helper function: ! fill vec with N elements from parms, starting at position ii !========================================================================== @@ -19,7 +12,7 @@ SUBROUTINE daisie_fill1d (vec, DIMP, parms, II) II = II + 1 vec(I) = parms(II) ENDDO - + END SUBROUTINE daisie_fill1d !========================================================================== @@ -31,10 +24,10 @@ MODULE daisie_dimmod ! length of the vector - decided in R-code INTEGER :: N INTEGER :: kk - + ! 1 parameter vectors with unknown length - DOUBLE PRECISION, ALLOCATABLE :: P(:) - + DOUBLE PRECISION, ALLOCATABLE :: P(:) + ! Boolean: will become TRUE if the parameters have a value LOGICAL :: initialised = .FALSE. @@ -48,39 +41,39 @@ END MODULE daisie_dimmod !========================================================================== SUBROUTINE daisie_initmod (steadyparms) - USE daisie_dimmod + USE daisie_dimmod IMPLICIT NONE EXTERNAL steadyparms INTEGER, PARAMETER :: nparsmall = 2 ! constant-length parameters - + DOUBLE PRECISION parms(nparsmall) - COMMON /XCBPar/parms ! common block + COMMON /XCBPar/parms ! common block ! Set the fixed parameters obtained from R CALL steadyparms(nparsmall, parms) -! first parameter has the length of the vector +! first parameter has the length of the vector N = INT(parms(1) + 1e-6) kk = INT(parms(2) + 1e-6) ! Allocate variable size arrays (state variables, derivatives and parameters) - IF (ALLOCATED(P)) DEALLOCATE(P) + IF (ALLOCATED(P)) DEALLOCATE(P) ALLOCATE(P(5 * (N + 4 + 2 * kk))) initialised = .FALSE. - + END SUBROUTINE daisie_initmod - + !========================================================================== !========================================================================== ! Dynamic routine: name of this function as passed by "func" argument ! variable parameter values are passed via yout !========================================================================== !========================================================================== - + SUBROUTINE daisie_runmod (neq, t, Conc, dConc, yout, ip) USE daisie_dimmod IMPLICIT NONE @@ -95,7 +88,6 @@ SUBROUTINE daisie_runmod (neq, t, Conc, dConc, yout, ip) DOUBLE PRECISION :: laavec(N + 4 + 2 * kk),lacvec(N + 4 + 2 * kk) DOUBLE PRECISION :: muvec(N + 4 + 2 * kk),gamvec(N + 4 + 2 * kk) DOUBLE PRECISION :: nn(N + 4 + 2 * kk) - DOUBLE PRECISION :: FF1, FFF ! parameters - named here DOUBLE PRECISION rn(2) @@ -108,7 +100,7 @@ SUBROUTINE daisie_runmod (neq, t, Conc, dConc, yout, ip) IF (.NOT. Initialised) THEN ! check memory allocated to output variables - IF (ip(1) < 1) CALL rexit("nout not large enough") + IF (ip(1) < 1) CALL rexit("nout not large enough") ! save parameter values in yout ii = ip(1) ! Start of parameter values @@ -117,7 +109,7 @@ SUBROUTINE daisie_runmod (neq, t, Conc, dConc, yout, ip) ENDIF ! dynamics - + ! xx1 = c(0,0,x[1:lx],0) ! xx2 = c(0,0,x[(lx + 1):(2 * lx)],0) ! xx3 = x[2 * lx + 1] @@ -165,55 +157,223 @@ SUBROUTINE daisie_runmod (neq, t, Conc, dConc, yout, ip) nn(I) = P(I + 4 * (N + 4 + 2 * kk)) ENDDO -! 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] +! dx3 <- 0 + + DO I = 1, N + dConc(I) = & + laavec(il1(I) + 1) * xx2(ix1(I)) + & + lacvec(il4(I) + 1) * xx2(ix4(I)) + & + muvec(il2(I) + 1) * xx2(ix3(I)) + & + lacvec(il1(I)) * nn(in1(I)) * xx1(ix1(I)) + & + muvec(il2(I)) * nn(in2ix2(I)) * xx1(in2ix2(I)) - & + (muvec(il3in3(I)) + lacvec(il3in3(I))) * & + nn(il3in3(I)) * xx1(ix3(I)) - & + gamvec(il3in3(I)) * xx1(ix3(I)) + dConc(N + I) = & + gamvec(il3in3(I)) * xx1(ix3(I)) + & + lacvec(il1(I) + 1) * nn(in1(I)) * xx2(ix1(I)) + & + muvec(il2(I) + 1) * nn(in2ix2(I)) * xx2(in2ix2(I)) - & + (muvec(il3in3(I) + 1) + lacvec(il3in3(I) + 1)) * & + nn(il3in3(I) + 1) * xx2(ix3(I)) - & + laavec(il3in3(I) + 1) * xx2(ix3(I)) + dConc(2 * N + 1) = 0 - DO I = 1, N - FF1 = laavec(il1(I) + 1) * xx2(ix1(I)) - FF1 = FF1 + lacvec(il4(I) + 1) * xx2(ix4(I)) - FF1 = FF1 + muvec(il2(I) + 1) * xx2(ix3(I)) - FF1 = FF1 + lacvec(il1(I)) * nn(in1(I)) * xx1(ix1(I)) - FF1 = FF1 + muvec(il2(I)) * nn(in2ix2(I)) * xx1(in2ix2(I)) - FFF = (muvec(il3in3(I)) + lacvec(il3in3(I))) - FF1 = FF1 - FFF * nn(il3in3(I)) * xx1(ix3(I)) - FF1 = FF1 - gamvec(il3in3(I)) * xx1(ix3(I)) - dConc(I) = FF1 - FF1 = gamvec(il3in3(I)) * xx1(ix3(I)) - FF1 = FF1 + lacvec(il1(I) + 1) * nn(in1(I)) * xx2(ix1(I)) - FF1 = FF1 + muvec(il2(I)+1) * nn(in2ix2(I)) * xx2(in2ix2(I)) - FFF = muvec(il3in3(I) + 1) + lacvec(il3in3(I) + 1) - FF1 = FF1 - FFF * nn(il3in3(I) + 1) * xx2(ix3(I)) - FF1 = FF1 - laavec(il3in3(I) + 1) * xx2(ix3(I)) - dConc(N + I) = FF1 ENDDO - IF(kk .EQ. 1) THEN - dConc(1) = dConc(1) + laavec(il3in3(1)) * xx3 - dConc(2) = dConc(2) + 2 * lacvec(il3in3(1)) * xx3 + + END SUBROUTINE daisie_runmod + +!========================================================================== +!========================================================================== +! Dynamic routine: name of this function as passed by "func" argument +! variable parameter values are passed via yout +!========================================================================== +!========================================================================== + + SUBROUTINE daisie_runmod1 (neq, t, Conc, dConc, yout, ip) + USE daisie_dimmod + IMPLICIT NONE + +!......................... declaration section............................. + INTEGER :: neq, ip(*), i, ii + DOUBLE PRECISION :: t, Conc(4 * N), dConc(4 * N), yout(*) + DOUBLE PRECISION :: xx1(N + 3), xx2(N + 3), xx3(N + 3), xx4(N + 3) + INTEGER :: il1(N), il2(N), il3in3(N), il4(N) + INTEGER :: in1(N), in2ix2(N), in4ix1(N) + INTEGER :: ix3(N), ix4(N) + DOUBLE PRECISION :: laavec(N + 4 + 2 * kk),lacvec(N + 4 + 2 * kk) + DOUBLE PRECISION :: muvec(N + 4 + 2 * kk),gamvec(N + 4 + 2 * kk) + DOUBLE PRECISION :: nn(N + 4 + 2 * kk) + DOUBLE PRECISION :: FF1, FFF + +! parameters - named here + DOUBLE PRECISION rn(2) + COMMON /XCBPar/rn + +! local variables + CHARACTER(len=100) msg + +!............................ statements .................................. + + IF (.NOT. Initialised) THEN + ! check memory allocated to output variables + IF (ip(1) < 1) CALL rexit("nout not large enough") + + ! save parameter values in yout + ii = ip(1) ! Start of parameter values + CALL daisie_fill1d(P, 5 * (N + 4 + 2 * kk), yout, ii) ! ii is updated in fill1d + Initialised = .TRUE. ! to prevent from initialising more than once ENDIF -! dx3 = -(laavec[il3in3[1]] + lacvec[il3in3[1]] + gamvec[il3in3[1]] + muvec[il3in3[1]]) * xx3 +! dynamics + +! xx1 <- c(0,0,x[1:lx],0) +! xx2 <- c(0,0,x[(lx + 1):(2 * lx)],0) +! xx3 <- c(0,0,x[(2 * lx + 1):(3 * lx)],0) +! xx4 <- c(0,0,x[(3 * lx + 1):(4 * lx)],0) +! nil2lx = 3:(lx + 2) +! il1 = nil2lx+kk-1 +! il2 = nil2lx+kk+1 +! il3 = nil2lx+kk +! il4 = nil2lx+kk-2 +! in1 = nil2lx+2*kk-1 +! in2 = nil2lx+1 +! in3 = nil2lx+kk +! in4 = nil2lx-1 +! ix1 = nil2lx-1 +! ix2 = nil2lx+1 +! ix3 = nil2lx +! ix4 = nil2lx-2 +! new: il3in3 = il3 = in3 +! new: in2ix2 = in2 = ix2 +! new: in4ix1 = in4 = ix1 + + xx1(1) = 0 + xx1(2) = 0 + xx2(1) = 0 + xx2(2) = 0 + xx3(1) = 0 + xx3(2) = 0 + xx4(1) = 0 + xx4(2) = 0 + DO I = 3, N + 2 + xx1(I) = Conc(I - 2) + xx2(I) = Conc(N + I - 2) + xx3(I) = Conc(2 * N + I - 2) + xx4(I) = Conc(3 * N + I - 2) + il1(I - 2) = I + kk - 1 + il2(I - 2) = I + kk + 1 + il3in3(I - 2) = I + kk + il4(I - 2) = I + kk - 2 + in1(I - 2) = I + 2 * kk - 1 + in2ix2(I - 2) = I + 1 + in4ix1(I - 2) = I - 1 + ix3(I - 2) = I + ix4(I - 2) = I - 2 + ENDDO + xx1(N + 3) = 0 + xx2(N + 3) = 0 + xx3(N + 3) = 0 + xx4(N + 3) = 0 + + DO I = 1, N + 4 + 2 * kk + laavec(I) = P(I) + lacvec(I) = P(I + N + 4 + 2 * kk) + muvec(I) = P(I + 2 * (N + 4 + 2 * kk)) + gamvec(I) = P(I + 3 * (N + 4 + 2 * kk)) + nn(I) = P(I + 4 * (N + 4 + 2 * kk)) + ENDDO + +! dx1 <- lacvec[il1] * xx1[ix1] + +! laavec[il1 + 1] * xx2[ix1] + +! lacvec[il4 + 1] * xx2[ix4] + +! muvec[il2] * nn[in2] * xx1[ix2] + +! muvec[il3 + 1] * xx2[ix3] + +! -(muvec[il3] + lacvec[il3]) * nn[in3] * xx1[ix3] + +! -gamvec[il3] * xx1[ix3] + + DO I = 1, N + dConc(I) = & + lacvec(il1(I)) * nn(in1(I)) * xx1(in4ix1(I)) + & + laavec(il1(I) + 1) * xx2(in4ix1(I)) + & + lacvec(il4(I) + 1) * xx2(ix4(I)) + & + muvec(il2(I)) * nn(in2ix2(I)) * xx1(in2ix2(I)) + & + muvec(il3in3(I) + 1) * xx2(ix3(I)) - & + (muvec(il3in3(I)) + lacvec(il3in3(I))) * & + nn(il3in3(I)) * xx1(ix3(I)) - & + gamvec(il3in3(I)) * xx1(ix3(I)) + +! dx2 <- gamvec[il3] * xx1[ix3] + +! gamvec[il3] * xx3[ix3] + +! gamvec[il3 + 1] * xx4[ix3] + +! 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] + + dConc(N + I) = & + gamvec(il3in3(I)) * xx1(ix3(I)) + & + gamvec(il3in3(I)) * xx3(ix3(I)) + & + gamvec(il3in3(I) + 1) * xx4(ix3(I)) + & + lacvec(il1(I) + 1) * nn(in1(I)) * xx2(in4ix1(I)) + & + muvec(il2(I) + 1) * nn(in2ix2(I)) * xx2(in2ix2(I)) - & + (muvec(il3in3(I) + 1) + lacvec(il3in3(I) + 1)) * & + nn(il3in3(I) + 1) * xx2(ix3(I)) - & + laavec(il3in3(I) + 1) * xx2(ix3(I)) + +! dx3 <- lacvec[il1] * nn[in1] * xx3[ix1] + +! laavec[il1 + 1] * xx4[ix1] + +! lacvec[il4 + 1] * xx4[ix4] + +! muvec[il2] * nn[in2] * xx3[ix2] + +! muvec[il3 + 1] * xx4[ix3] + +! -(lacvec[il3] + muvec[il3]) * nn[in3] * xx3[ix3] + +! -gamvec[il3] * xx3[ix3] + + dConc(2 * N + I) = & + lacvec(il1(I)) * nn(in1(I)) * xx3(in4ix1(I)) + & + laavec(il1(I) + 1) * xx4(in4ix1(I)) + & + lacvec(il4(I) + 1) * xx4(ix4(I)) + & + muvec(il2(I)) * nn(in2ix2(I)) * xx3(in2ix2(I)) + & + muvec(il3in3(I) + 1) * xx4(ix3(I)) - & + (lacvec(il3in3(I)) + muvec(il3in3(I))) * & + nn(il3in3(I)) * xx3(ix3(I)) - & + gamvec(il3in3(I)) * xx3(ix3(I)) + +! dx4 <- lacvec[il1 + 1] * nn[in1] * xx4[ix1] + +! muvec[il2 + 1] * nn[in2] * xx4[ix2] + +! -(lacvec[il3 + 1] + muvec[il3 + 1]) * nn[in3 + 1] * xx4[ix3] + +! -gamvec[il3 + 1] * xx4[ix3] + + dConc(3 * N + I) = & + lacvec(il1(I) + 1) * nn(in1(I)) * xx4(in4ix1(I)) + & + muvec(il2(I) + 1) * nn(in2ix2(I)) * xx4(in2ix2(I)) - & + (lacvec(il3in3(I) + 1) + muvec(il3in3(I) + 1)) * & + nn(il3in3(I) + 1) * xx4(ix3(I)) - & + gamvec(il3in3(I) + 1) * xx4(ix3(I)) + + ENDDO + + END SUBROUTINE daisie_runmod1 - FFF = laavec(il3in3(1)) + lacvec(il3in3(1)) - FFF = FFF + gamvec(il3in3(1)) + muvec(il3in3(1)) - dConc(2 * N + 1) = -1 * FFF * xx3 - END SUBROUTINE daisie_runmod - !========================================================================== !========================================================================== ! Dynamic routine: name of this function as passed by "func" argument ! variable parameter values are passed via yout !========================================================================== !========================================================================== - + SUBROUTINE daisie_runmod2 (neq, t, Conc, dConc, yout, ip) USE daisie_dimmod IMPLICIT NONE @@ -228,7 +388,6 @@ SUBROUTINE daisie_runmod2 (neq, t, Conc, dConc, yout, ip) DOUBLE PRECISION :: laavec(N + 4 + 2 * kk),lacvec(N + 4 + 2 * kk) DOUBLE PRECISION :: muvec(N + 4 + 2 * kk),gamvec(N + 4 + 2 * kk) DOUBLE PRECISION :: nn(N + 4 + 2 * kk) - DOUBLE PRECISION :: FF1, FFF ! parameters - named here DOUBLE PRECISION rn(2) @@ -241,7 +400,7 @@ SUBROUTINE daisie_runmod2 (neq, t, Conc, dConc, yout, ip) IF (.NOT. Initialised) THEN ! check memory allocated to output variables - IF (ip(1) < 1) CALL rexit("nout not large enough") + IF (ip(1) < 1) CALL rexit("nout not large enough") ! save parameter values in yout ii = ip(1) ! Start of parameter values @@ -250,7 +409,7 @@ SUBROUTINE daisie_runmod2 (neq, t, Conc, dConc, yout, ip) ENDIF ! dynamics - + ! xx1 = c(0,0,x[1:lx],0) ! xx2 = c(0,0,x[(lx + 1):(2 * lx)],0) ! xx3 = c(0,0,x[(2 * lx + 1):(3 * lx)],0) @@ -303,47 +462,56 @@ SUBROUTINE daisie_runmod2 (neq, t, Conc, dConc, yout, ip) nn(I) = P(I + 4 * (N + 4 + 2 * kk)) ENDDO -! dx1 = (laavec[il3] * xx3[ix3] + 2 * lacvec[il1] * xx3[in4ix1]) * (kk == 1) + -! laavec[il1 + 1] * xx2[in4ix1] + lacvec[il4 + 1] * xx2[ix4] + muvec[il2 + 1] * xx2[ix3] -! lacvec[il1] * nn[in1] * xx1[in4ix1] + muvec[il2] * nn[in2] * xx1[in2ix2] + -! -(muvec[il3] + lacvec[il3]) * nn[in3] * xx1[ix3] - gamvec[il3] * xx1[ix3] -! dx2 = gamvec[il3] * xx1[ix3] + -! lacvec[il1 + 1] * nn[in1] * xx2[in4ix1] + muvec[il2 + 1] * nn[in2] * xx2[in2ix2] + +! 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] + +! dx2 <- gamvec[il3] * xx1[ix3] + +! 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] -! dx3 = lacvec[il1] * nn[in4ix1] * xx3[in4ix1] + muvec[il2] * nn[in2] * xx3[in2ix2] + + +! 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] - - DO I = 1, N - FF1 = laavec(il1(I) + 1) * xx2(in4ix1(I)) - FF1 = FF1 + lacvec(il4(I) + 1) * xx2(ix4(I)) - FF1 = FF1 + muvec(il2(I) + 1) * xx2(ix3(I)) - FF1 = FF1 + lacvec(il1(I)) * nn(in1(I)) * xx1(in4ix1(I)) - FF1 = FF1 + muvec(il2(I)) * nn(in2ix2(I)) * xx1(in2ix2(I)) - FFF = muvec(il3in3(I)) + lacvec(il3in3(I)) - FF1 = FF1 - FFF * nn(il3in3(I)) * xx1(ix3(I)) - FF1 = FF1 - gamvec(il3in3(I)) * xx1(ix3(I)) - FFF = 0 + + DO I = 1, N + dConc(I) = 0 IF(kk .EQ. 1) THEN - FFF = laavec(il3in3(I)) * xx3(ix3(I)) - FFF = FFF + 2 * lacvec(il1(I)) * xx3(in4ix1(I)) + dConc(I) = laavec(il3in3(I)) * xx3(ix3(I)) + & + 2 * lacvec(il1(I)) * xx3(in4ix1(I)) ENDIF - dConc(I) = FF1 + FFF - FF1 = gamvec(il3in3(I)) * xx1(ix3(I)) - FF1 = FF1 + lacvec(il1(I) + 1) * nn(in1(I)) * xx2(in4ix1(I)) - FF1 = FF1 + muvec(il2(I)+1) * nn(in2ix2(I)) * xx2(in2ix2(I)) - FFF = muvec(il3in3(I) + 1) + lacvec(il3in3(I) + 1) - FF1 = FF1 - FFF * nn(il3in3(I) + 1) * xx2(ix3(I)) - FF1 = FF1 - laavec(il3in3(I) + 1) * xx2(ix3(I)) - dConc(N + I) = FF1 - FF1 = lacvec(il1(I)) * nn(in4ix1(I)) * xx3(in4ix1(I)) - FF1 = FF1 + muvec(il2(I)) * nn(in2ix2(I)) * xx3(in2ix2(I)) - FFF = lacvec(il3in3(I)) + muvec(il3in3(I)) - FF1 = FF1 - FFF * nn(il3in3(I)) * xx3(ix3(I)) - FF1 = FF1-(laavec(il3in3(I))+gamvec(il3in3(I)))*xx3(ix3(I)) - dConc(2 * N + I) = FF1 + dConc(I) = dConc(I) + & + laavec(il1(I) + 1) * xx2(in4ix1(I)) + & + lacvec(il4(I) + 1) * xx2(ix4(I)) + & + muvec(il2(I) + 1) * xx2(ix3(I)) + & + lacvec(il1(I)) * nn(in1(I)) * xx1(in4ix1(I)) + & + muvec(il2(I)) * nn(in2ix2(I)) * xx1(in2ix2(I)) - & + (muvec(il3in3(I)) + lacvec(il3in3(I))) * & + nn(il3in3(I)) * xx1(ix3(I)) - & + gamvec(il3in3(I)) * xx1(ix3(I)) + dConc(N + I) = & + gamvec(il3in3(I)) * xx1(ix3(I)) + & + lacvec(il1(I) + 1) * nn(in1(I)) * xx2(in4ix1(I)) + & + muvec(il2(I) + 1) * nn(in2ix2(I)) * xx2(in2ix2(I)) - & + (muvec(il3in3(I) + 1) + lacvec(il3in3(I) + 1)) * & + nn(il3in3(I) + 1) * xx2(ix3(I)) - & + laavec(il3in3(I) + 1) * xx2(ix3(I)) + dConc(2 * N + I) = & + lacvec(il1(I)) * nn(in4ix1(I)) * xx3(in4ix1(I)) + & + muvec(il2(I)) * nn(in2ix2(I)) * xx3(in2ix2(I)) - & + (lacvec(il3in3(I)) + muvec(il3in3(I))) * & + nn(il3in3(I)) * xx3(ix3(I)) - & + (laavec(il3in3(I)) + gamvec(il3in3(I))) * xx3(ix3(I)) ENDDO END SUBROUTINE daisie_runmod2 - + diff --git a/src/DAISIE_odeint.h b/src/DAISIE_odeint.h index 905cdaf6..4bea02d3 100644 --- a/src/DAISIE_odeint.h +++ b/src/DAISIE_odeint.h @@ -6,31 +6,119 @@ #define STRICT_R_HEADERS +#include + +// [[Rcpp::plugins(cpp14)]] +// [[Rcpp::depends(BH)]] + +#include + +// Provide Forward Declarations +namespace Rcpp { + + namespace traits{ + + // Setup non-intrusive extension via template specialization for + // 'ublas' class boost::numeric::ublas + + // Support for wrap + template SEXP wrap(const boost::numeric::ublas::vector & obj); + + // Support for as + template class Exporter< boost::numeric::ublas::vector >; + + } +} + + #include #include #include +#include #include +// boost::numeric::ublas wrapping from: +// https://gallery.rcpp.org/articles/custom-templated-wrap-and-as-for-seamingless-interfaces/ +namespace Rcpp { + + namespace traits{ + + // Defined wrap case + template SEXP wrap(const boost::numeric::ublas::vector & obj){ + const int RTYPE = Rcpp::traits::r_sexptype_traits::rtype ; + + return Rcpp::Vector< RTYPE >(obj.begin(), obj.end()); + }; + + + // Defined as< > case + template class Exporter< boost::numeric::ublas::vector > { + typedef typename boost::numeric::ublas::vector OUT ; + + // Convert the type to a valid rtype. + const static int RTYPE = Rcpp::traits::r_sexptype_traits< T >::rtype ; + Rcpp::Vector vec; + + public: + Exporter(SEXP x) : vec(x) { + if (TYPEOF(x) != RTYPE) + throw std::invalid_argument("Wrong R type for mapped 1D array"); + } + OUT get() { + // Need to figure out a way to perhaps do a pointer pass? + OUT x(vec.size()); + std::copy(vec.begin(), vec.end(), x.begin()); // have to copy data + return x; + } + }; + + } + +} + + using namespace Rcpp; using namespace boost::numeric::odeint; // type of the ode state -using state_type = std::vector; +using state_type = boost::numeric::ublas::vector; + namespace daisie_odeint { - template - inline void do_integrate(double atol, double rtol, Rhs rhs, state_type& y, double t0, double t1) - { - integrate_adaptive(make_controlled(atol, rtol), rhs, y, t0, t1, 0.1 * (t1 - t0)); - } +extern double abm_factor; - // wrapper around odeint::integrate +template +inline void do_integrate(double atol, double rtol, Rhs rhs, state_type& y, double t0, double t1) +{ + integrate_adaptive(make_controlled(atol, rtol), rhs, y, t0, t1, 0.1 * (t1 - t0)); +} + + +template +inline void abm(Rhs rhs, state_type& y, double t0, double t1) +{ + auto abm = adams_bashforth_moulton{}; + abm.initialize(rhs, y, t0, abm_factor * (t1 - t0)); + integrate_const(abm, rhs, y, t0, t1, abm_factor * (t1 - t0)); +} + + +template +inline void ab(Rhs rhs, state_type& y, double t0, double t1) +{ + auto ab = adams_bashforth{}; + ab.initialize(rhs, y, t0, abm_factor * (t1 - t0)); + integrate_const(ab, rhs, y, t0, t1, abm_factor * (t1 - t0)); +} + + +// wrapper around odeint::integrate // maps runtime stepper name -> compiletime odeint::stepper type template inline void integrate( @@ -56,6 +144,34 @@ namespace daisie_odeint { using stepper_t = bulirsch_stoer; integrate_adaptive(stepper_t(atol, rtol), rhs, y, t0, t1, 0.1 * (t1 - t0)); } + else if (0 == stepper.compare(0, stepper.size() - 2, "odeint::adams_bashforth")) { + const char steps = stepper.back(); + switch (steps) { + case '1': ab<1>(rhs, y, t0, t1); break; + case '2': ab<2>(rhs, y, t0, t1); break; + case '3': ab<3>(rhs, y, t0, t1); break; + case '4': ab<4>(rhs, y, t0, t1); break; + case '5': ab<5>(rhs, y, t0, t1); break; + case '6': ab<6>(rhs, y, t0, t1); break; + case '7': ab<7>(rhs, y, t0, t1); break; + case '8': ab<8>(rhs, y, t0, t1); break; + default: throw std::runtime_error("DAISIE_odeint_helper::integrate: unsupported steps for admam_bashforth_moulton"); + } + } + else if (0 == stepper.compare(0, stepper.size() - 2, "odeint::adams_bashforth_moulton")) { + const char steps = stepper.back(); + switch (steps) { + case '1': abm<1>(rhs, y, t0, t1); break; + case '2': abm<2>(rhs, y, t0, t1); break; + case '3': abm<3>(rhs, y, t0, t1); break; + case '4': abm<4>(rhs, y, t0, t1); break; + case '5': abm<5>(rhs, y, t0, t1); break; + case '6': abm<6>(rhs, y, t0, t1); break; + case '7': abm<7>(rhs, y, t0, t1); break; + case '8': abm<8>(rhs, y, t0, t1); break; + default: throw std::runtime_error("DAISIE_odeint_helper::integrate: unsupported steps for admam_bashforth_moulton"); + } + } else { throw std::runtime_error("DAISIE_odeint_helper::integrate: unknown stepper"); } diff --git a/src/R_init_DAISIE.c b/src/R_init_DAISIE.c index 63a18e76..7dd3fcf5 100644 --- a/src/R_init_DAISIE.c +++ b/src/R_init_DAISIE.c @@ -7,12 +7,14 @@ extern void F77_NAME(daisie_fill1d)(double *vec, int *DIMP, double *parms, int *II); extern void F77_NAME(daisie_initmod)(void (*steadyparms)(int *, double *)); extern void F77_NAME(daisie_runmod)(int *neq, double *t, double *Conc, double *dConc, double *yout, int *ip); +extern void F77_NAME(daisie_runmod1)(int *neq, double *t, double *Conc, double *dConc, double *yout, int *ip); extern void F77_NAME(daisie_runmod2)(int *neq, double *t, double *Conc, double *dConc, double *yout, int *ip); static const R_FortranMethodDef FortranEntries[] = { {"daisie_fill1d", (DL_FUNC) &F77_NAME(daisie_fill1d), 4}, {"daisie_initmod", (DL_FUNC) &F77_NAME(daisie_initmod), 1}, {"daisie_runmod", (DL_FUNC) &F77_NAME(daisie_runmod), 6}, + {"daisie_runmod1", (DL_FUNC) &F77_NAME(daisie_runmod1), 6}, {"daisie_runmod2", (DL_FUNC) &F77_NAME(daisie_runmod2), 6}, {NULL, NULL, 0} }; diff --git a/tests/testthat/test-DAISIE_ML1.R b/tests/testthat/test-DAISIE_ML1.R index 9540c333..2a11295f 100644 --- a/tests/testthat/test-DAISIE_ML1.R +++ b/tests/testthat/test-DAISIE_ML1.R @@ -30,7 +30,7 @@ test_that("use", { K = 906.6501180193454, gamma = 0.0173458887696076, lambda_a = 3.677789527566334, - loglik = -64.22005531779574, + loglik = -64.2199684450019, df = 5L, conv = 0L ) diff --git a/tests/testthat/test-DAISIE_ML2.R b/tests/testthat/test-DAISIE_ML2.R index 733a55a3..3d33840a 100644 --- a/tests/testthat/test-DAISIE_ML2.R +++ b/tests/testthat/test-DAISIE_ML2.R @@ -1,6 +1,6 @@ test_that("use", { - skip_if(Sys.getenv("CI") == "", message = "Run only on CI") + utils::data(Macaronesia_datalist, package = "DAISIE") invisible(capture.output( tested_MLE <- DAISIE_ML2( @@ -27,29 +27,29 @@ test_that("use", { )) expected_MLE <- data.frame( lambda_c = c(0.0, - 0.333180238362478, + 0.315015693879739, 0.0, 0.0), - mu = c(2.718309102127384, - 2.718309102127384, - 2.718309102127384, - 2.718309102127384), + mu = c(2.67928476648249, + 2.67928476648249, + 2.67928476648249, + 2.67928476648249), K = c(Inf, Inf, Inf, Inf), - 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), + gamma = c(0.141748213992604, + 0.338044376412124, + 0.141748213992604, + 0.141748213992604), + lambda_a = c(1.10809210249609, + 1.10809210249609, + 1.10809210249609 , + 1.10809210249609), + loglik = c(-409.60117813893, + -409.60117813893, + -409.60117813893, + -409.60117813893), 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 deleted file mode 100644 index d721bd88..00000000 --- a/tests/testthat/test-DAISIE_ML3.R +++ /dev/null @@ -1,85 +0,0 @@ -test_that("use", { - skip_if(Sys.getenv("CI") == "", message = "Run only on CI") - # 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) - 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] - ) - 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, - 0.133766934, - 0.0, - 0.0), - mu = c(1.0531518319999997, - 1.0531518319999997, - 1.0531518319999997, - 1.0531518319999997), - K = c(Inf, - Inf, - Inf, - Inf), - gamma = c(0.052148978999999998, - 0.152763178999999999, - 0.052148978999999998, - 0.052148978999999998), - lambda_a = c(0.51293901099999994, - 0.51293901099999994, - 0.51293901099999994, - 0.51293901099999994), - loglik = c(-454.93478332906614, - -454.93478332906614, - -454.93478332906614, - -454.93478332906614), - df = c(5L, 5L, 5L, 5L), - conv = c(0L, 0L, 0L, 0L) - ) - -}) - -test_that("abuse", { - skip_if(Sys.getenv("CI") == "", message = "Run only on CI") - expect_error(tested_MLE <- DAISIE:::DAISIE_ML2( - datalist = "nonsense", - 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) - )) -}) diff --git a/tests/testthat/test-DAISIE_MW_ML.R b/tests/testthat/test-DAISIE_MW_ML.R index 1b78311b..875dfc93 100644 --- a/tests/testthat/test-DAISIE_MW_ML.R +++ b/tests/testthat/test-DAISIE_MW_ML.R @@ -2,27 +2,12 @@ test_that("DAISIE_MW_ML produces correct output", { skip_if(Sys.getenv("CI") == "", message = "Run only on CI") - utils::data(archipelagos41) - M19_Nature_parameters <- c( - 4.007380e-02, - 0.000000e+00, - 1.945657e+00, - 1.504297e-01, - Inf, - 0.000000e+00, - 67.2564367200001, - 2.936351e-01, - 5.909687e-02, - 3.826885e-01, - 2.651078e-02, - -3651.6624531664, - 8.000000e+00, - 0.000000e+00 - ) + + utils::data(archipelagos41) invisible(capture.output( - M19_computation <- DAISIE::DAISIE_MW_ML( + M19_tested <- DAISIE::DAISIE_MW_ML( datalist = archipelagos41, initparsopt = c( 0.040073803, @@ -48,10 +33,27 @@ test_that("DAISIE_MW_ML produces correct output", { ) )) + M19_Nature_expected <- c( + 0.040073803, + 0.0, + 1.945656546, + 0.150429656, + Inf, + 0.0, + 67.2564367200001, + 0.293635061, + 0.059096872, + 0.382688527, + 0.026510781, + -3651.48307905, + 8, + 0 + ) + + testthat::expect_equal( - M19_Nature_parameters, - as.numeric(M19_computation), - tol = 0.000001 + M19_Nature_expected, + as.numeric(M19_tested) ) }) diff --git a/tests/testthat/test-DAISIE_format_IW.R b/tests/testthat/test-DAISIE_format_IW.R index 634899d7..b2506e96 100644 --- a/tests/testthat/test-DAISIE_format_IW.R +++ b/tests/testthat/test-DAISIE_format_IW.R @@ -539,22 +539,27 @@ test_that("silent when species with two trait states with stt_all[2, ] <- c(0, 0, 1, 2, 0, 0, 0) 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) + brts_table[1, ] <- c(5.00000000000000, 0, 0, NA, NA) + brts_table[2, ] <- c(3.10634202528338, 1, 1, 1, NA) + brts_table[3, ] <- c(1.52330128016821, 2, 1, 1, NA) + brts_table[4, ] <- c(1.28012784155125, 2, 2, 1, NA) expected_IW_format[[1]][[1]] <- list(island_age = 5, not_present = 13, stt_all = stt_all, brts_table = brts_table) - expected_IW_format[[1]][[2]] <- list(branching_times = c(5.000000000, - 3.102613675), - stac = 2, - missing_species = 0) - expected_IW_format[[1]][[3]] <- list(branching_times = c(5.000000000, - 1.505629998, - 1.262456559), - stac = 2, - missing_species = 0) - expect_true(all.equal(formatted_IW_sim, expected_IW_format, tolerance = 1e-7)) + + expected_IW_format[[1]][[2]] <- list( + branching_times = c(5.00000000000000, + 3.10634202528338), + stac = 2, + missing_species = 0 + ) + expected_IW_format[[1]][[3]] <- list( + branching_times = c(5.00000000000000, + 1.52330128016821, + 1.28012784155125), + stac = 2, + missing_species = 0 + ) + expect_equal(formatted_IW_sim, expected_IW_format, tolerance = 1e-7) }) diff --git a/tests/testthat/test-DAISIE_loglik_CS.R b/tests/testthat/test-DAISIE_loglik_CS.R index 3b8aee1a..071c5036 100644 --- a/tests/testthat/test-DAISIE_loglik_CS.R +++ b/tests/testthat/test-DAISIE_loglik_CS.R @@ -45,6 +45,7 @@ test_that("DAISIE_loglik_CS_choice produces correct output for relaxed-rate test_that("DAISIE_loglik_CS_choice produces same output for CS_version = 0 (with M = 1) and CS_version = 1 ", { + skip_if(Sys.getenv("CI") == "", message = "Run only on CI") pars1 <- c(2.000, 2.700, 20.000, 0.009, 1.010) pars2 <- c(100, 11, 0, 0, NA, 0.0e+00, 1.0e-04, @@ -70,7 +71,7 @@ test_that("DAISIE_loglik_CS_choice produces same output for CS_version = 0 missnumspec = missnumspec, CS_version = CS_version) - expect_equal(loglik0,loglik1) + expect_equal(loglik0, loglik1) }) test_that("DAISIE_loglik_all produces correct output for relaxed-rate model", { @@ -108,8 +109,11 @@ test_that("DAISIE_loglik produces correct output", { testthat::expect_equal(output, -0.00347317077256095) }) -test_that("DAISIE_loglik_all produces same output for CS_version 0 and 1 with and without conditioning", { + +test_that("DAISIE_loglik_all produces same output for CS_version 0 and 1 with + and without conditioning", { skip_if(Sys.getenv("CI") == "", message = "Run only on CI") + utils::data(Galapagos_datalist) Galapagos_datalist2 <- Galapagos_datalist for(i in 2:9) { @@ -150,7 +154,7 @@ test_that("DAISIE_loglik_all produces same output for CS_version 0 and 1 with an CS_version = 1, abstolint = 1e-16, reltolint = 1e-10) - testthat::expect_equal(loglik_CS01, loglik_CS11, tol = 5E-6) + expect_equal(loglik_CS01, loglik_CS11, tol = 5E-6) }) test_that("DAISIE_loglik_CS_choice produces equivalent output for ODEINT RKCK54 @@ -163,25 +167,26 @@ test_that("DAISIE_loglik_CS_choice produces equivalent output for ODEINT RKCK54 0.0527, 0.0327, 0.0221, 0.1180, 0.0756, 0.0525, 0.0322, 0.0118) stac <- 2 missnumspec <- 0 - loglik1 <- DAISIE:::DAISIE_loglik_CS_choice(pars1 = pars1, - pars2 = pars2, - brts = brts, - stac = stac, - missnumspec = missnumspec) - loglik2 <- DAISIE:::DAISIE_loglik_CS_choice(pars1 = pars1, - pars2 = pars2, - brts = brts, - stac = stac, - missnumspec = missnumspec, - methode = "odeint::runge_kutta_cash_karp54") + loglik1 <- DAISIE:::DAISIE_loglik_CS_choice( + pars1 = pars1, + pars2 = pars2, + brts = brts, + stac = stac, + missnumspec = missnumspec + ) + loglik2 <- DAISIE:::DAISIE_loglik_CS_choice( + pars1 = pars1, + pars2 = pars2, + brts = brts, + stac = stac, + missnumspec = missnumspec, + methode = "odeint::runge_kutta_cash_karp54" + ) expect_equal(expected = loglik1, object = loglik2) }) - - -test_that("DAISIE_loglik_CS_choice produces equivalent output for ODEINT - RKF78 and deSolve lsodes", { - skip_if(Sys.getenv("CI") == "", message = "Run only on CI") +test_that("DAISIE_loglik_CS_choice produces equivalent + output for ODEINT RKF78 and deSolve lsodes", { 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) @@ -191,19 +196,29 @@ test_that("DAISIE_loglik_CS_choice produces equivalent output for ODEINT missnumspec <- 0 CS_version <- 0 # deSolve lsodes - loglik1 <- DAISIE_loglik_CS_choice(pars1 = pars1, - pars2 = pars2, - brts = brts, - stac = stac, - missnumspec = missnumspec, - CS_version = CS_version) + + loglik1 <- expect_silent( + DAISIE_loglik_CS_choice( + pars1 = pars1, + pars2 = pars2, + brts = brts, + stac = stac, + missnumspec = missnumspec, + CS_version = CS_version + ) + ) # odeint RKF78 - loglik2 <- DAISIE_loglik_CS_choice(pars1 = pars1, - pars2 = pars2, - brts = brts, - stac = stac, - missnumspec = missnumspec, - CS_version = CS_version, - methode = "odeint::runge_kutta_fehlberg78") + loglik2 <- expect_silent( + DAISIE_loglik_CS_choice( + pars1 = pars1, + pars2 = pars2, + brts = brts, + stac = stac, + missnumspec = missnumspec, + CS_version = CS_version, + methode = "odeint::runge_kutta_fehlberg78" + ) + ) expect_equal(expected = loglik1, object = loglik2) }) + diff --git a/tests/testthat/test-DAISIE_loglik_CS_M1.R b/tests/testthat/test-DAISIE_loglik_CS_M1.R new file mode 100644 index 00000000..f5a28301 --- /dev/null +++ b/tests/testthat/test-DAISIE_loglik_CS_M1.R @@ -0,0 +1,91 @@ +test_that("DAISIE_loglik_CS_M1 produces correct output",{ + + dataset <- list( + list(island_age = 1, not_present = 90), + list(branching_times = c(1.0, 0.99999), stac = 1, missing_species = 0), + list(branching_times = c(1.0, 0.999989), stac = 1, missing_species = 0), + list(branching_times = c(1.0, 0.5), stac = 1, missing_species = 0), + list(branching_times = c(1.0, 0.99999), stac = 5, missing_species = 0), + list(branching_times = c(1.0, 0.99989), stac = 5, missing_species = 0), + list(branching_times = c(1.0, 0.5), stac = 5, missing_species = 0), + list(branching_times = c(1.0, 0.99999, 0.000001), stac = 8, + missing_species = 0), + list(branching_times = c(1.0, 0.99989, 0.000001), stac = 8, + missing_species = 0), + list(branching_times = c(1.0, 0.5, 0.000001), stac = 8, + missing_species = 0), + list(branching_times = c(1.0, 0.99999, 0.000001), stac = 9, + missing_species = 0), + list(branching_times = c(1.0, 0.99989, 0.000001), stac = 9, + missing_species = 0), + list(branching_times = c(1.0, 0.5, 0.000001), stac = 9, + missing_species = 0) + ) + out_1 <- c() + out_2 <- c() + for (i in 2:length(dataset)) { + + invisible(capture.output(out_1[i] <- DAISIE_loglik_CS_M1( + brts = dataset[[i]]$branching_times, + stac = dataset[[i]]$stac, + missnumspec = dataset[[i]]$missing_species, + pars1 = c(0.1, 0.1, 10.0, 0.1, 0.1), + pars2 = c(1.0e+02, 1.1e+01, 0.0e+00, 1.0e+00, NA, 0.0e+00, 1.0e-04, + 1.0e-05, + 1.0e-07, 3.0e+03, 9.5e-01, 9.8e-01), + verbose = FALSE + ))) + invisible(capture.output(out_2[i] <- DAISIE_loglik_CS_M1( + brts = dataset[[i]]$branching_times, + stac = dataset[[i]]$stac, + missnumspec = dataset[[i]]$missing_species, + pars1 = c(0.5, 0.1, 10.0, 0.1, 0.1), + pars2 = c(1.0e+02, 1.1e+01, 0.0e+00, 1.0e+00, NA, 0.0e+00, 1.0e-04, + 1.0e-05, + 1.0e-07, 3.0e+03, 9.5e-01, 9.8e-01), + verbose = FALSE + ))) + } + out_1 <- out_1[-1] + out_2 <- out_2[-1] + expected_out_1 <- c(-2.49433738044967, + -2.49434644748344, + -3.09980594295824, + -5.46225301144597, + -5.46245593750427, + -6.77188944406840, + -2.49433857201992, + -2.49442924792230, + -3.09980813006604, + -5.46225301144714, + -5.46245593750544, + -6.77188944407274 + ) + expected_out_2 <- c( + -2.66327113437861, + -2.66327862430581, + -3.19372583697100, + -5.51403419701533, + -5.51422850883500, + -6.79274094715768, + -2.66327252801121, + -2.66334743269888, + -3.19372820911879, + -5.51403419701655, + -5.51422850883622, + -6.79274094716205 + ) + expect_equal(out_1, expected_out_1) + expect_equal(out_2, expected_out_2) + + # Max_ages at island age should be very close to max ages at very close to + # island age + expect_lt(out_1[1] - out_1[2], 1e-3) + expect_lt(out_1[4] - out_1[5], 1e-3) + expect_lt(out_1[7] - out_1[8], 1e-3) + expect_lt(out_1[10] - out_1[11], 1e-3) + expect_lt(out_2[1] - out_2[2], 1e-3) + expect_lt(out_2[4] - out_2[5], 1e-3) + expect_lt(out_2[7] - out_2[8], 1e-3) + expect_lt(out_2[10] - out_2[11], 1e-3) +}) diff --git a/tests/testthat/test-DAISIE_sim_core_cr.R b/tests/testthat/test-DAISIE_sim_core_cr.R index 0d5a8175..1a50caec 100644 --- a/tests/testthat/test-DAISIE_sim_core_cr.R +++ b/tests/testthat/test-DAISIE_sim_core_cr.R @@ -86,7 +86,7 @@ test_that("DAISIE_sim_core output is correct", { expect_true(is.numeric(sim_core$taxon_list[[1]]$branching_times)) expect_true(is.numeric(sim_core$taxon_list[[1]]$stac)) expect_true(is.numeric(sim_core$taxon_list[[1]]$missing_species)) - expect_true(length(sim_core$taxon_list) == 2) + expect_true(length(sim_core$taxon_list) == 5) expect_true("branching_times" %in% names(sim_core$taxon_list[[1]])) expect_true("stac" %in% names(sim_core$taxon_list[[1]])) expect_true("missing_species" %in% names(sim_core$taxon_list[[1]])) diff --git a/tests/testthat/test-integration_DAISIE.R b/tests/testthat/test-integration_DAISIE.R index d7978733..5c3fdaec 100644 --- a/tests/testthat/test-integration_DAISIE.R +++ b/tests/testthat/test-integration_DAISIE.R @@ -3,7 +3,7 @@ test_that("loglik Galapagos works", { rm(Galapagos_datalist) Galapagos_datalist_2types <- NULL rm(Galapagos_datalist_2types) - utils::data(Galapagos_datalist_2types, package = "DAISIE") + data(Galapagos_datalist_2types, package = "DAISIE") pars1 <- c( 0.195442017, 0.087959583, @@ -25,7 +25,7 @@ test_that("loglik Galapagos works", { test_that("loglik macaronesia 2 type works", { Macaronesia_datalist <- NULL rm(Macaronesia_datalist) - utils::data(Macaronesia_datalist, package = "DAISIE") + data(Macaronesia_datalist, package = "DAISIE") background <- c(0, 1.053151832, Inf, 0.052148979, 0.512939011) Canaries <- c(0.133766934, 1.053151832, Inf, 0.152763179, 0.512939011) pars1 <- rbind(background, Canaries, background, background) @@ -33,15 +33,15 @@ test_that("loglik macaronesia 2 type works", { loglik <- 0 for (i in seq_along(Macaronesia_datalist)) { loglik <- loglik + DAISIE_loglik_all(pars1[i, ], - pars2, - Macaronesia_datalist[[i]], - methode = "lsodes") + pars2, + Macaronesia_datalist[[i]], + methode = "lsodes") } - testthat::expect_equal(loglik, -452.1156003854809) + expect_equal(loglik, -449.921430187808) }) test_that("clade specific rate-shift loglik works", { - utils::data(Galapagos_datalist, package = "DAISIE") + data(Galapagos_datalist, package = "DAISIE") pars1 <- c(0.2, 0.1, Inf, 0.001, 0.3, 0.2, 0.1, Inf, 0.001, 0.3, 1) pars2 <- c(40, 11, 0, 0) SR_loglik_CS <- DAISIE_SR_loglik_CS( @@ -57,13 +57,13 @@ test_that("clade specific rate-shift loglik works", { datalist = Galapagos_datalist, methode = "ode45", CS_version = 1) - testthat::expect_equal(SR_loglik_CS, loglik_CS) + expect_equal(SR_loglik_CS, loglik_CS) }) test_that("IW and CS loglik is same when K = Inf", { skip_if(Sys.getenv("CI") == "" && !(Sys.getenv("USERNAME") == "rampa"), message = "Run only on CI") - utils::data(Galapagos_datalist, package = "DAISIE") + data(Galapagos_datalist, package = "DAISIE") pars1 <- c(0.35, 0.3, Inf, 0.001, 0.3) pars2 <- c(80, 11, 1, 0) Galapagos_datalist_IW <- Galapagos_datalist @@ -71,34 +71,64 @@ test_that("IW and CS loglik is same when K = Inf", { 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 <- add_brt_table(Galapagos_datalist_IW) - invisible(capture.output( - loglik_IW <- DAISIE_loglik_IW( - pars1 = pars1, - pars2 = pars2, - datalist = Galapagos_datalist_IW, - methode = "ode45") - )) - invisible(capture.output( - loglik_IW2 <- DAISIE_loglik_IW( + + Galapagos_datalist_IW <- DAISIE:::add_brt_table(Galapagos_datalist_IW) + loglik_IW <- DAISIE_loglik_IW( pars1 = pars1, pars2 = pars2, datalist = Galapagos_datalist_IW, - methode = "odeint::runge_kutta_fehlberg78") - )) + methode = "ode45" + ) + + loglik_IW2 <- DAISIE_loglik_IW( + pars1 = pars1, + pars2 = pars2, + datalist = Galapagos_datalist_IW, + methode = "odeint::runge_kutta_fehlberg78" + ) + + + loglik_CS <- DAISIE_loglik_CS( + pars1 = pars1, + pars2 = pars2, + datalist = Galapagos_datalist_IW, + methode = "ode45", + CS_version = 1 + ) + + expect_equal(loglik_IW, loglik_IW2, tol = 5E-6) + expect_equal(loglik_IW, loglik_CS, tol = 5E-6) +}) + +test_that("DAISIE_ML simple case works", { + skip_if(Sys.getenv("CI") == "" && !(Sys.getenv("USERNAME") == "rampa"), + message = "Run only on CI") + expected_mle <- data.frame( + 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) + invisible(capture.output( - loglik_CS <- DAISIE_loglik_CS( - pars1 = pars1, - pars2 = pars2, - datalist = Galapagos_datalist_IW, - methode = "ode45", - CS_version = 1) + tested_mle <- 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(loglik_IW, loglik_IW2, tol = 5E-6) - testthat::expect_equal(loglik_IW, loglik_CS, tol = 5E-6) + expect_equal(expected_mle, tested_mle) }) -testthat::test_that("DAISIE_ML simple case works", { +test_that("DAISIE_ML simple case works with zero probability of initial presence", { skip_if(Sys.getenv("CI") == "" && !(Sys.getenv("USERNAME") == "rampa"), message = "Run only on CI") expected_mle <- data.frame( @@ -107,32 +137,118 @@ testthat::test_that("DAISIE_ML simple case works", { K = 2992.207701921788, gamma = 0.00937711049761019, lambda_a = 0.9993246958280274, + prob_init_pres = 0, loglik = -75.99266304738612, df = 5L, conv = 0L ) utils::data(Galapagos_datalist) - cat("\n") + + invisible(capture.output( + tested_mle <- DAISIE_ML( + datalist = Galapagos_datalist, + initparsopt = c(2.5, 2.7, 20, 0.009, 1.01), + ddmodel = 11, + idparsopt = 1:5, + parsfix = 0, + idparsfix = 6 + ) + )) + expect_equal(expected_mle, tested_mle) +}) + +test_that("DAISIE_ML simple case works with nonzero probability of initial presence", { + skip_if(Sys.getenv("CI") == "" && !(Sys.getenv("USERNAME") == "rampa"), + message = "Run only on CI") + expected_mle <- data.frame( + lambda_c = 2.58373135630384, + mu = 2.70882802751483, + K = 2992.20770192179, + gamma = 0.00937711049761019, + lambda_a = 0.999324695828027, + prob_init_pres = 0.1, + loglik = -75.9926628720867, + df = 5L, + conv = 0L + ) + utils::data(Galapagos_datalist) + invisible(capture.output( tested_mle <- DAISIE_ML( datalist = Galapagos_datalist, initparsopt = c(2.5, 2.7, 20, 0.009, 1.01), ddmodel = 11, idparsopt = 1:5, + parsfix = 0.1, + idparsfix = 6 + ) + )) + expect_equal(expected_mle, tested_mle) +}) + +test_that("DAISIE_ML simple case works with estimating probability of initial + presence", { + skip_if(Sys.getenv("CI") == "" && !(Sys.getenv("USERNAME") == "rampa"), + message = "Run only on CI") + if (identical(Sys.getenv("OS"), "Windows_NT")) { + expected_mle <- data.frame( + lambda_c = 2.53430497145461, + mu = 2.66658569091753, + K = 2136343.97554965, + gamma = 0.00930345848936764, + lambda_a = 1.0119011474385, + prob_init_pres = 3.21939792431987e-10, + loglik = -75.9925548510873, + df = 6L, + conv = 0L + ) + } else if (identical(Sys.info()["sysname"], c(sysname = "Darwin"))) { + expected_mle <- data.frame( + lambda_c = 2.5330538395353, + mu = 2.66544727106831, + K = 8477857.16185865, + gamma = 0.00929919368081989, + lambda_a = 1.01226940726093, + prob_init_pres = 2.08612789348566e-08, + loglik = -75.9925781743204, + df = 6L, + conv = 0L + ) + } else { + expected_mle <- data.frame( + lambda_c = 2.53432108511347, + mu = 2.66677757261811, + K = 2155153.9420102, + gamma = 0.00930305175196706, + lambda_a = 1.01184784588089, + prob_init_pres = 2.33936540081158e-10, + loglik = -75.9925542005831, + df = 6L, + conv = 0L + ) + } + + utils::data(Galapagos_datalist) + invisible(capture.output( + tested_mle <- DAISIE_ML( + datalist = Galapagos_datalist, + initparsopt = c(2.5, 2.7, 20, 0.009, 1.01, 0.001), + ddmodel = 11, + idparsopt = 1:6, parsfix = NULL, idparsfix = NULL ) )) - testthat::expect_equal(expected_mle, tested_mle) + expect_equal(tested_mle, expected_mle) }) test_that("The parameter choice for 2type DAISIE_ML works", { Galapagos_datalist_2types <- NULL rm(Galapagos_datalist_2types) - utils::data(Galapagos_datalist_2types, package = "DAISIE") + data(Galapagos_datalist_2types, package = "DAISIE") set.seed(1) # MLE and high tolerance for speed-up - cat("\n") + invisible(capture.output( fit <- DAISIE_ML( datalist = Galapagos_datalist_2types, @@ -146,26 +262,55 @@ test_that("The parameter choice for 2type DAISIE_ML works", { maxiter = 30 ) )) - testthat::expect_equal(fit$loglik, -74.7557, tol = 1E-3) + expect_equal(fit$loglik, -74.7557, tol = 1E-3) }) test_that("conditioning works", { # Cond 0 ## 1 type - utils::data(Galapagos_datalist, package = "DAISIE") + data(Galapagos_datalist, package = "DAISIE") pars1_1type_cond0 <- c(0.2, 0.1, Inf, 0.001, 0.3) + #pars1_1type_cond0 <- c(0, 0, Inf, 1, 0) pars2_1type_cond0 <- c(40, 11, 0, 0) - loglik_CS_1type_cond0 <- DAISIE_loglik_CS( + res1 <- DAISIE_loglik_CS( pars1 = pars1_1type_cond0, pars2 = pars2_1type_cond0, datalist = Galapagos_datalist, methode = "ode45", CS_version = 1 ) + res2 <- DAISIE_loglik_CS( + pars1 = pars1_1type_cond0, + pars2 = pars2_1type_cond0, + datalist = Galapagos_datalist, + methode = "deSolve_R::ode45", + CS_version = 1 + ) + res3 <- loglik_CS_1type_cond0 <- DAISIE_loglik_CS( + pars1 = pars1_1type_cond0, + pars2 = pars2_1type_cond0, + datalist = Galapagos_datalist, + methode = "odeint::runge_kutta_fehlberg78", + CS_version = 1 + ) + + testthat::expect_equal(res1, res3) + testthat::expect_equal(res2, res3, tol = 1E-4) testthat::expect_equal(loglik_CS_1type_cond0, -96.49069330275196) + # Status of colonist: 0, Parameters: 0.200000 0.100000 Inf 0.001000 0.300000 , Loglikelihood: -0.003424 + # Status of colonist: 1, Parameters: 0.200000 0.100000 Inf 0.001000 0.300000 , Loglikelihood: -6.494398 + # Status of colonist: 4, Parameters: 0.200000 0.100000 Inf 0.001000 0.300000 , Loglikelihood: -7.113751 + # Status of colonist: 2, Parameters: 0.200000 0.100000 Inf 0.001000 0.300000 , Loglikelihood: -31.251817 + # Status of colonist: 2, Parameters: 0.200000 0.100000 Inf 0.001000 0.300000 , Loglikelihood: -14.421388 + # Status of colonist: 2, Parameters: 0.200000 0.100000 Inf 0.001000 0.300000 , Loglikelihood: -8.594293 + # Status of colonist: 2, Parameters: 0.200000 0.100000 Inf 0.001000 0.300000 , Loglikelihood: -10.599996 + # Status of colonist: 1, Parameters: 0.200000 0.100000 Inf 0.001000 0.300000 , Loglikelihood: -6.494398 + # Status of colonist: 2, Parameters: 0.200000 0.100000 Inf 0.001000 0.300000 , Loglikelihood: -8.123768 + # Parameters: 0.200000 0.100000 Inf 0.001000 0.300000 , Loglikelihood: -96.490693 + ## 2 type - utils::data(Galapagos_datalist_2types, package = "DAISIE") + data(Galapagos_datalist_2types, package = "DAISIE") pars1_2type_cond0 <- c( 0.195442017, 0.087959583, @@ -185,11 +330,11 @@ test_that("conditioning works", { pars2_2type_cond0, Galapagos_datalist_2types ) - testthat::expect_equal(loglik_CS_2type_cond0, -61.70281911731144) + expect_equal(loglik_CS_2type_cond0, -61.7028188767349) # Cond 1 ## 1 type - utils::data(Galapagos_datalist, package = "DAISIE") + data(Galapagos_datalist, package = "DAISIE") pars1_1type_cond1 <- c(0.2, 0.1, Inf, 0.001, 0.3) pars2_1type_cond1 <- c(40, 11, 1, 0) loglik_CS_1type_cond1 <- DAISIE_loglik_CS( @@ -199,10 +344,10 @@ test_that("conditioning works", { methode = 'ode45', CS_version = 1 ) - testthat::expect_equal(loglik_CS_1type_cond1, -96.45757823017264) + expect_equal(loglik_CS_1type_cond1, -96.45757823017264) ## 2 type - utils::data(Galapagos_datalist_2types, package = "DAISIE") + data(Galapagos_datalist_2types, package = "DAISIE") pars1_2type_cond1 <- c( 0.195442017, 0.087959583, @@ -222,7 +367,7 @@ test_that("conditioning works", { pars2_2type_cond1, Galapagos_datalist_2types ) - testthat::expect_equal(loglik_CS_2type_cond1, -61.4375956792401) + expect_equal(loglik_CS_2type_cond1, -61.4375954386635) # Cond 5 ## 1 type @@ -236,10 +381,10 @@ test_that("conditioning works", { methode = 'ode45', CS_version = 1 ) - testthat::expect_equal(loglik_CS_1type_cond5, -95.14000237210625) + expect_equal(loglik_CS_1type_cond5, -95.14000237210625) ## 2 type - utils::data(Galapagos_datalist_2types, package = "DAISIE") + data(Galapagos_datalist_2types, package = "DAISIE") pars1_2type_cond5 <- c( 0.195442017, 0.087959583, @@ -259,5 +404,5 @@ test_that("conditioning works", { pars2_2type_cond5, Galapagos_datalist_2types ) - testthat::expect_equal(loglik_CS_2type_cond5, -61.3735196464293) + expect_equal(loglik_CS_2type_cond5, -61.3735194058527) }) diff --git a/tests/testthat/test-update_rates.R b/tests/testthat/test-update_rates.R index 82e577cf..502c243a 100644 --- a/tests/testthat/test-update_rates.R +++ b/tests/testthat/test-update_rates.R @@ -57,11 +57,12 @@ test_that("update area-dependent rates is silent and gives correct output", { max_area = 1.0, current_area = 0.5, proportional_peak_t = 0.5, - total_island_age = 1.0, + total_island_age = 1.2, sea_level_amplitude = 0, sea_level_frequency = 0, island_gradient_angle = 0) hyper_pars <- create_hyper_pars(d = 0.2, x = 0.1) + peak <- calc_peak(total_time = 1, area_pars = area_pars) expect_silent(rates <- DAISIE:::update_rates( timeval = 0, total_time = 1, @@ -77,7 +78,8 @@ test_that("update area-dependent rates is silent and gives correct output", { K = K, num_spec = 0, num_immigrants = 0, - mainland_n = 1)) + mainland_n = 1, + peak = peak)) expect_true(are_rates(rates)) expected_rates <- list(immig_rate = 0, ext_rate = 0,