diff --git a/.github/workflows/test-coverage.yaml b/.github/workflows/test-coverage.yaml index 952a4f7..fdb6800 100644 --- a/.github/workflows/test-coverage.yaml +++ b/.github/workflows/test-coverage.yaml @@ -1,50 +1,63 @@ -# Workflow derived from https://github.com/r-lib/actions/tree/v2/examples -# Need help debugging build failures? Start at https://github.com/r-lib/actions#where-to-find-help on: push: + branches: [master, develop] pull_request: + branches: [master, develop] -name: test-coverage +name: test-coverage-new + +permissions: read-all jobs: test-coverage: runs-on: ubuntu-latest - if: "contains(github.event.head_commit.message, '[run ci]') || (github.ref == 'refs/heads/master' || github.ref == 'refs/heads/develop')" - env: GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} steps: - - uses: actions/checkout@v3 + - uses: actions/checkout@v4 - uses: r-lib/actions/setup-r@v2 with: use-public-rspm: true - + + - name: Install devtools + run: install.packages("devtools") + shell: Rscript {0} + - uses: r-lib/actions/setup-r-dependencies@v2 with: - extra-packages: any::covr + extra-packages: any::covr, any::xml2 needs: coverage - name: Test coverage run: | - covr::codecov( + cov <- covr::package_coverage( quiet = FALSE, clean = FALSE, - install_path = file.path(Sys.getenv("RUNNER_TEMP"), "package") + install_path = file.path(normalizePath(Sys.getenv("RUNNER_TEMP"), winslash = "/"), "package") ) + covr::to_cobertura(cov) shell: Rscript {0} + - uses: codecov/codecov-action@v5 + with: + fail_ci_if_error: ${{ github.event_name != 'pull_request' && true || false }} + files: ./cobertura.xml + plugin: noop + disable_search: true + token: ${{ secrets.CODECOV_TOKEN }} + - name: Show testthat output if: always() run: | ## -------------------------------------------------------------------- - find ${{ runner.temp }}/package -name 'testthat.Rout*' -exec cat '{}' \; || true + find '${{ runner.temp }}/package' -name 'testthat.Rout*' -exec cat '{}' \; || true shell: bash - name: Upload test results if: failure() - uses: actions/upload-artifact@v3 + uses: actions/upload-artifact@v4 with: name: coverage-test-failures - path: ${{ runner.temp }}/package + path: ${{ runner.temp }}/package \ No newline at end of file diff --git a/DESCRIPTION b/DESCRIPTION index 0a22769..c6686ca 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,8 +1,7 @@ Package: DDD Type: Package Title: Diversity-Dependent Diversification -Version: 5.2.3 -Date: 2024-11-26 +Version: 5.2.4 Depends: R (>= 3.5.0) Imports: deSolve, diff --git a/R/RcppExports.R b/R/RcppExports.R index 7497dee..4087da9 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -8,6 +8,13 @@ dd_integrate_bw_odeint <- function(ry, times, pars, atol, rtol, stepper) { #' @useDynLib DDD NULL +dd_integrate_log_odeint <- function(ry, times, pars, atol, rtol, stepper) { + .Call('_DDD_dd_integrate_log_odeint', PACKAGE = 'DDD', ry, times, pars, atol, rtol, stepper) +} + +#' @useDynLib DDD +NULL + dd_integrate_odeint <- function(ry, times, pars, atol, rtol, stepper) { .Call('_DDD_dd_integrate_odeint', PACKAGE = 'DDD', ry, times, pars, atol, rtol, stepper) } diff --git a/R/bd_ML.R b/R/bd_ML.R index c069309..136f238 100644 --- a/R/bd_ML.R +++ b/R/bd_ML.R @@ -57,7 +57,7 @@ #' @param changeloglikifnoconv if TRUE the loglik will be set to -Inf if ML #' does not converge #' @param optimmethod Method used in optimization of the likelihood. Current -#' default is 'subplex'. Alternative is 'simplex' (default of previous +#' default is 'simplex'. Alternative is 'subplex' (default of previous #' versions) #' @param num_cycles the number of cycles of opimization. If set at Inf, it will #' do as many cycles as needed to meet the tolerance set for the target function. @@ -102,7 +102,7 @@ bd_ML = function(brts, tol = c(1E-3, 1E-4, 1E-6), maxiter = 1000 * round((1.25)^length(idparsopt)), changeloglikifnoconv = FALSE, - optimmethod = 'subplex', + optimmethod = 'simplex', num_cycles = 1, methode = 'odeint::runge_kutta_cash_karp54', verbose = FALSE) diff --git a/R/dd_KI_ML.R b/R/dd_KI_ML.R index b03b2c3..19d02da 100644 --- a/R/dd_KI_ML.R +++ b/R/dd_KI_ML.R @@ -78,7 +78,7 @@ #' @param changeloglikifnoconv if TRUE the loglik will be set to -Inf if ML #' does not converge #' @param optimmethod Method used in optimization of the likelihood. Current -#' default is 'subplex'. Alternative is 'simplex' (default of previous +#' default is 'simplex'. Alternative is 'subplex' (default of previous #' versions) #' @param num_cycles the number of cycles of opimization. If set at Inf, it will #' do as many cycles as needed to meet the tolerance set for the target function. @@ -140,7 +140,7 @@ dd_KI_ML = function(brtsM, tol = c(1E-3, 1E-4, 1E-6), maxiter = 1000 * round((1.25)^length(idparsopt)), changeloglikifnoconv = FALSE, - optimmethod = 'subplex', + optimmethod = 'simplex', num_cycles = 1, methode = 'analytical', correction = TRUE, diff --git a/R/dd_LR.R b/R/dd_LR.R index 9d96042..7244180 100644 --- a/R/dd_LR.R +++ b/R/dd_LR.R @@ -60,7 +60,7 @@ #' @param changeloglikifnoconv if TRUE the loglik will be set to -Inf if ML #' does not converge #' @param optimmethod Method used in optimization of the likelihood. Current -#' default is 'subplex'. Alternative is 'simplex' (default of previous +#' default is 'simplex'. Alternative is 'subplex' (default of previous #' versions) #' @param methode The method used to solve the master equation, default is #' 'analytical' which uses matrix exponentiation; alternatively numerical ODE @@ -120,7 +120,7 @@ dd_LR = function( tol = c(1E-3,1E-4,1E-6), maxiter = 2000, changeloglikifnoconv = FALSE, - optimmethod = 'subplex', + optimmethod = 'simplex', methode = 'analytical' ) { diff --git a/R/dd_ML.R b/R/dd_ML.R index e6ebd05..7c1fd4c 100644 --- a/R/dd_ML.R +++ b/R/dd_ML.R @@ -91,15 +91,20 @@ parsfixdefault = function(ddmodel,brts,missnumspec,idparsopt) #' @param btorph Sets whether the likelihood is for the branching times (0) or #' the phylogeny (1) #' @param soc Sets whether stem or crown age should be used (1 or 2) -#' @param tol Sets the tolerances in the optimization. Consists of: \cr reltolx -#' = relative tolerance of parameter values in optimization \cr reltolf = -#' relative tolerance of function value in optimization \cr abstolx = absolute -#' tolerance of parameter values in optimization +#' @param tol Sets the tolerances in the optimization. Consists of:\cr +#' reltolx = relative tolerance of parameter values in optimization \cr +#' reltolf = relative tolerance of function value in optimization \cr +#' abstolx = absolute tolerance of parameter values in optimization +#' @param tolint Sets the tolerance of the numerical integration. Consists of: \cr +#' absoltint = absolute tolerance and \cr +#' reltolint = relative tolerance. +#' @param probs_threshold Sets the threshold of the probability below which +#' logarithmic integration must be used. Default is 0. #' @param maxiter Sets the maximum number of iterations in the optimization #' @param changeloglikifnoconv if TRUE the loglik will be set to -Inf if ML #' does not converge #' @param optimmethod Method used in optimization of the likelihood. Current -#' default is 'subplex'. Alternative is 'simplex' (default of previous +#' default is 'simplex'. Alternative is 'subplex' (default of previous #' versions) #' @param num_cycles the number of cycles of opimization. If set at Inf, it will #' do as many cycles as needed to meet the tolerance set for the target function. @@ -145,9 +150,11 @@ dd_ML = function( btorph = 1, soc = 2, tol = c(1E-3, 1E-4, 1E-6), + tolint = c(1E-10,1E-8), + probs_threshold = 0, maxiter = 1000 * round((1.25)^length(idparsopt)), changeloglikifnoconv = FALSE, - optimmethod = 'subplex', + optimmethod = 'simplex', num_cycles = 1, methode = 'analytical', verbose = FALSE) @@ -183,7 +190,7 @@ dd_ML = function( trparsopt[which(initparsopt == Inf)] = 1 trparsfix = parsfix/(1 + parsfix) trparsfix[which(parsfix == Inf)] = 1 - pars2 = c(res,ddmodel,cond,btorph,verbose,soc,tol,maxiter) + pars2 = c(res,ddmodel,cond,btorph,verbose,soc,tol,maxiter,abstolint = tolint[1],reltolint = tolint[2],probs_threshold = probs_threshold) optimpars = c(tol,maxiter) initloglik = dd_loglik_choosepar(trparsopt = trparsopt,trparsfix = trparsfix,idparsopt = idparsopt,idparsfix = idparsfix,pars2 = pars2,brts = brts,missnumspec = missnumspec, methode = methode) cat("The loglikelihood for the initial parameter values is",initloglik,"\n") diff --git a/R/dd_MS_ML.R b/R/dd_MS_ML.R index 8be2fd8..c7ad62d 100644 --- a/R/dd_MS_ML.R +++ b/R/dd_MS_ML.R @@ -78,7 +78,7 @@ #' @param changeloglikifnoconv if TRUE the loglik will be set to -Inf if ML #' does not converge #' @param optimmethod Method used in optimization of the likelihood. Current -#' default is 'subplex'. Alternative is 'simplex' (default of previous +#' default is 'simplex'. Alternative is 'subplex' (default of previous #' versions) #' @param num_cycles the number of cycles of opimization. If set at Inf, it will #' do as many cycles as needed to meet the tolerance set for the target function. @@ -136,7 +136,7 @@ dd_MS_ML = function(brtsM, tol = c(1E-3, 1E-4, 1E-6), maxiter = 1000 * round((1.25)^length(idparsopt)), changeloglikifnoconv = FALSE, - optimmethod = 'subplex', + optimmethod = 'simplex', num_cycles = 1, methode = 'ode45', correction = FALSE, diff --git a/R/dd_SR_ML.R b/R/dd_SR_ML.R index f9d8893..4059c42 100644 --- a/R/dd_SR_ML.R +++ b/R/dd_SR_ML.R @@ -71,7 +71,7 @@ #' @param changeloglikifnoconv if TRUE the loglik will be set to -Inf if ML #' does not converge #' @param optimmethod Method used in optimization of the likelihood. Current -#' default is 'subplex'. Alternative is 'simplex' (default of previous +#' default is 'simplex'. Alternative is 'subplex' (default of previous #' versions) #' @param num_cycles the number of cycles of opimization. If set at Inf, it will #' do as many cycles as needed to meet the tolerance set for the target function. @@ -128,7 +128,7 @@ dd_SR_ML = function(brts, tol = c(1E-3, 1E-4, 1E-6), maxiter = 1000 * round((1.25)^length(idparsopt)), changeloglikifnoconv = FALSE, - optimmethod = 'subplex', + optimmethod = 'simplex', num_cycles = 1, methode = 'analytical', verbose = FALSE) diff --git a/R/dd_loglik.R b/R/dd_loglik.R index 335fa77..3fe4318 100644 --- a/R/dd_loglik.R +++ b/R/dd_loglik.R @@ -109,6 +109,11 @@ dd_loglik_test = function(pars1,pars2,brts,missnumspec,methode = 'analytical',rh #' on screen (1) or not (0) #' \cr \cr \code{pars2[6]} sets whether the first data #' point is stem age (1) or crown age (2) +#' \cr \cr \code{pars2[7]} sets the absolute tolerance of the integration +#' \cr \cr \code{pars2[8]} sets the relative tolerance of the integration +#' \cr \cr \code{pars2[9]} sets threshold on the number of species for which +#' logarithmic integration will be used. If the last three elements are not set +#' ddefault values will be used of 1E-10, 1E-8 and Inf. #' @param brts A set of branching times of a phylogeny, all positive #' @param missnumspec The number of species that are in the clade but missing #' in the phylogeny @@ -132,21 +137,21 @@ dd_loglik_test = function(pars1,pars2,brts,missnumspec,methode = 'analytical',rh #' @export dd_loglik dd_loglik = function(pars1,pars2,brts,missnumspec,methode = 'analytical') { - if(pars2[3] == 3) - { - rhs_func_name = 'dd_loglik_bw_rhs' - } else - { - rhs_func_name = 'dd_loglik_rhs' - } - if(methode == 'analytical') - { - out = dd_loglik2(pars1,pars2,brts,missnumspec) - } else - { - out = dd_loglik1(pars1,pars2,brts,missnumspec,methode = methode,rhs_func_name = rhs_func_name) - } - return(out) + if(pars2[3] == 3) + { + rhs_func_name = 'dd_loglik_bw_rhs' + } else + { + rhs_func_name = 'dd_loglik_rhs' + } + if(methode == 'analytical') + { + out = dd_loglik2(pars1,pars2,brts,missnumspec) + } else + { + out = dd_loglik1(pars1,pars2,brts,missnumspec,methode = methode,rhs_func_name = rhs_func_name) + } + return(out) } dd_loglik1 = function(pars1,pars2,brts,missnumspec,methode = 'odeint::runge_kutta_cash_karp54',rhs_func_name = 'dd_loglik_rhs') @@ -181,8 +186,9 @@ dd_loglik1 = function(pars1,pars2,brts,missnumspec,methode = 'odeint::runge_kutt { loglik = bd_loglik(pars1[1:(2 + (K < Inf))],c(2*(mu == 0 & K < Inf),pars2[3:6]),brts,missnumspec) } else { - abstol = 1e-10 - reltol = 1e-8 + if(is.na(pars2['abstolint'])) abstolint <- 1e-10 else abstolint <- pars2['abstolint'] + if(is.na(pars2['reltolint'])) reltolint <- 1e-8 else reltolint <- pars2['reltolint'] + if(is.na(pars2['probs_threshold'])) probs_threshold <- 0 else probs_threshold <- pars2['probs_threshold'] brts = -sort(abs(as.numeric(brts)),decreasing = TRUE) if(sum(brts == 0) == 0) { @@ -209,10 +215,21 @@ dd_loglik1 = function(pars1,pars2,brts,missnumspec,methode = 'odeint::runge_kutt { probs = rep(0,lx) probs[1] = 1 # change if other species at stem/crown age + lx_old <- lx + logtrafo <- FALSE for(k in 2:(S + 2 - soc)) { k1 = k + (soc - 2) - y = dd_integrate(probs,brts[(k-1):k],rhs_func_name,c(pars1,k1,ddep),rtol = reltol,atol = abstol,method = methode) + probs <- probs[1:lx] + if(probs[1] > 0 & probs[1] < probs_threshold) { + lx <- min(lx,which(probs == 0) - 1) + probs <- probs[1:lx] + logtrafo <- TRUE + rhs_func_name <- 'dd_loglik_log_rhs' + probs <- log(probs) + } + + y = dd_integrate(probs,brts[(k-1):k],rhs_func_name,c(pars1,k1,ddep),rtol = reltolint,atol = abstolint,method = methode) probs = y[2,2:(lx+1)] if(is.na(sum(probs)) && pars1[2]/pars1[1] < 1E-4 && missnumspec == 0) { @@ -222,9 +239,18 @@ dd_loglik1 = function(pars1,pars2,brts,missnumspec,methode = 'odeint::runge_kutt } if(k < (S + 2 - soc)) { - probs = flavec(ddep,la,mu,K,r,lx,k1) * probs # speciation event + if(logtrafo) { + fac <- flavec(ddep,la,mu,K,r,lx,k1) + lx <- min(lx,which(fac == 0) - 1) + probs <- probs[1:lx] + fac <- fac[1:lx] + probs <- log(fac) + probs + } else + { + probs = flavec(ddep,la,mu,K,r,lx,k1) * probs # speciation event + cp <- check_probs(loglik,probs,verbose); loglik <- cp[[1]]; probs <- cp[[2]]; + } } - cp <- check_probs(loglik,probs,verbose); loglik <- cp[[1]]; probs <- cp[[2]]; } } else { probs = rep(0,lx + 1) @@ -232,7 +258,7 @@ dd_loglik1 = function(pars1,pars2,brts,missnumspec,methode = 'odeint::runge_kutt for(k in (S + 2 - soc):2) { k1 = k + (soc - 2) - y = dd_integrate(probs,-brts[k:(k-1)],rhs_func_name,c(pars1,k1,ddep),rtol = reltol,atol = abstol,method = methode) + y = dd_integrate(probs,-brts[k:(k-1)],rhs_func_name,c(pars1,k1,ddep),rtol = reltolint,atol = abstolint,method = methode) probs = y[2,2:(lx+2)] if(k > soc) { @@ -241,14 +267,21 @@ dd_loglik1 = function(pars1,pars2,brts,missnumspec,methode = 'odeint::runge_kutt cp <- check_probs(loglik,probs[1:lx],verbose); loglik <- cp[[1]]; probs[1:lx] <- cp[[2]]; } } - if(probs[1 + missnumspec] <= 0 | loglik == -Inf | is.na(loglik) | is.nan(loglik)) + if((!logtrafo & probs[1 + missnumspec] <= 0) | loglik == -Inf | is.na(loglik) | is.nan(loglik)) { if(verbose) cat('Probabilities smaller than 0 or other numerical problems are encountered in final result.\n') loglik = -Inf } else { - loglik = loglik + (cond != 3 | soc == 1) * log(probs[1 + (cond != 3) * missnumspec]) - lgamma(S + missnumspec + 1) + lgamma(S + 1) + lgamma(missnumspec + 1) + if(!logtrafo) { + loglik = loglik + (cond != 3 | soc == 1) * log(probs[1 + (cond != 3) * missnumspec]) - lgamma(S + missnumspec + 1) + lgamma(S + 1) + lgamma(missnumspec + 1) + } else + { + loglik <- loglik + (cond != 3 | soc == 1) * probs[1 + (cond != 3) * missnumspec] - lgamma(S + missnumspec + 1) + lgamma(S + 1) + lgamma(missnumspec + 1) + rhs_func_name <- 'dd_loglik_rhs' + } logliknorm = 0 + lx <- lx_old if(cond == 1 | cond == 2) { probsn = rep(0,lx) @@ -256,14 +289,16 @@ dd_loglik1 = function(pars1,pars2,brts,missnumspec,methode = 'odeint::runge_kutt k = soc t1 = brts[1] t2 = brts[S + 2 - soc] - y = dd_integrate(probsn,c(t1,t2),rhs_func_name,c(pars1,k,ddep),rtol = reltol,atol = abstol,method = methode); + y = dd_integrate(probsn,c(t1,t2),rhs_func_name,c(pars1,k,ddep),rtol = reltolint,atol = abstolint,method = methode); probsn = y[2,2:(lx+1)] if(soc == 1) { aux = 1:lx } if(soc == 2) { aux = (2:(lx+1)) * (3:(lx+2))/6 } probsc = probsn/aux - cp <- check_probs(logliknorm,probsc,verbose); logliknorm <- cp[[1]]; probsc <- cp[[2]]; - if(cond == 1) { logliknorm = logliknorm + log(sum(probsc)) } - if(cond == 2) { logliknorm = logliknorm + log(probsc[S + missnumspec - soc + 1])} + #cp <- check_probs(logliknorm,probsc,verbose); logliknorm <- cp[[1]]; probsc <- cp[[2]]; + #if(cond == 1) { logliknorm = logliknorm + log(sum(probsc)) } + #if(cond == 2) { logliknorm = logliknorm + log(probsc[S + missnumspec - soc + 1])} + if(cond == 1) { logliknorm = log(sum(probsc)) } + if(cond == 2) { logliknorm = log(probsc[S + missnumspec - soc + 1])} } if(cond == 3) { @@ -271,14 +306,14 @@ dd_loglik1 = function(pars1,pars2,brts,missnumspec,methode = 'odeint::runge_kutt probsn[S + missnumspec + 1] = 1 TT = 1e14 # max(1,1/abs(la - mu)) * 1E+10 * max(abs(brts)) # make this more efficient later - y = dd_integrate(probsn,c(0,TT),rhs_func_name,c(pars1,0,ddep),rtol = reltol,atol = abstol,method = methode) + y = dd_integrate(probsn,c(0,TT),rhs_func_name,c(pars1,0,ddep),rtol = reltolint,atol = abstolint,method = methode) logliknorm = log(y[2,lx + 2]) if(soc == 2) { probsn = rep(0,lx + 1) probsn[1:lx] = probs[1:lx] probsn = c(flavec(ddep,la,mu,K,r,lx,1),1) * probsn # speciation event - y = dd_integrate(probsn,c(max(abs(brts)),TT),rhs_func_name,c(pars1,1,ddep),rtol = reltol,atol = abstol,method = methode) + y = dd_integrate(probsn,c(max(abs(brts)),TT),rhs_func_name,c(pars1,1,ddep),rtol = reltolint,atol = abstolint,method = methode) logliknorm = logliknorm - log(y[2,lx + 2]) } } @@ -313,151 +348,151 @@ dd_loglik1 = function(pars1,pars2,brts,missnumspec,methode = 'odeint::runge_kutt dd_loglik2 = function(pars1,pars2,brts,missnumspec) { -if(length(pars2) == 4) -{ + if(length(pars2) == 4) + { pars2[5] = 0 pars2[6] = 2 -} -ddep = pars2[2] -cond = pars2[3] -btorph = pars2[4] -verbose <- pars2[5] -soc = pars2[6] -if(cond == 3) -{ + } + ddep = pars2[2] + cond = pars2[3] + btorph = pars2[4] + verbose <- pars2[5] + soc = pars2[6] + if(cond == 3) + { soc = 2 -} -la = pars1[1] -mu = pars1[2] -K = pars1[3] -if(ddep == 5) -{ + } + la = pars1[1] + mu = pars1[2] + K = pars1[3] + if(ddep == 5) + { r = pars1[4] -} else -{ + } else + { r = 0 -} -if(ddep == 1 | ddep == 5) -{ + } + if(ddep == 1 | ddep == 5) + { lx = min(max(1 + missnumspec,1 + ceiling(la/(la - mu) * (r + 1) * K)),ceiling(pars2[1])) -} else if(ddep == 1.3) -{ + } else if(ddep == 1.3) + { lx = min(ceiling(K),ceiling(pars2[1])) -} else { + } else { lx = round(pars2[1]) -} -if((ddep == 1) & ((mu == 0 & missnumspec == 0 & floor(K) != ceiling(K) & la > 0.05) | K == Inf)) -{ + } + if((ddep == 1) & ((mu == 0 & missnumspec == 0 & floor(K) != ceiling(K) & la > 0.05) | K == Inf)) + { loglik = bd_loglik(pars1[1:(2 + (K < Inf))],c(2*(mu == 0 & K < Inf),pars2[3:6]),brts,missnumspec) -} else { -abstol = 1e-16 -reltol = 1e-10 -brts = -sort(abs(as.numeric(brts)),decreasing = TRUE) -if(sum(brts == 0) == 0) -{ - brts[length(brts) + 1] = 0 -} -S = length(brts) + (soc - 2) -if(min(pars1) < 0) -{ - loglik = -Inf -} else { -if((mu == 0 & (ddep == 2 | ddep == 2.1 | ddep == 2.2)) | (la == 0 & (ddep == 4 | ddep == 4.1 | ddep == 4.2)) | (la <= mu)) -{ - if(verbose) cat("These parameter values cannot satisfy lambda(N) = mu(N) for a positive and finite N.\n") - loglik = -Inf -} else { - if(((ddep == 1 | ddep == 5) & ceiling(la/(la - mu) * (r + 1) * K) < (S + missnumspec)) | ((ddep == 1.3) & ((S + missnumspec) > ceiling(K)))) + } else { + if(is.na(pars2['abstolint'])) abstolint <- 1e-16 else abstolint <- pars2['abstolint'] + if(is.na(pars2['reltolint'])) reltolint <- 1e-10 else reltolint <- pars2['reltolint'] + brts = -sort(abs(as.numeric(brts)),decreasing = TRUE) + if(sum(brts == 0) == 0) + { + brts[length(brts) + 1] = 0 + } + S = length(brts) + (soc - 2) + if(min(pars1) < 0) { - loglik = -Inf + loglik = -Inf } else { - loglik = (btorph == 0) * lgamma(S) - if(cond != 3) - { - probs = rep(0,lx) - probs[1] = 1 # change if other species at stem/crown age - for(k in 2:(S + 2 - soc)) - { - k1 = k + (soc - 2) - #y = deSolve::ode(probs,brts[(k-1):k],rhs_func,c(pars1,k1,ddep),rtol = reltol,atol = abstol,method = methode) - #probs2 = y[2,2:(lx+1)] - probs = dd_loglik_M(pars1,lx,k1,ddep,tt = abs(brts[k] - brts[k-1]),probs) - if(is.na(sum(probs)) && pars1[2]/pars1[1] < 1E-4 && missnumspec == 0) - { - loglik = dd_loglik_high_lambda(pars1 = pars1,pars2 = pars2,brts = brts) - if(verbose) cat('High lambda approximation has been applied.\n') - return(loglik) - } - if(k < (S + 2 - soc)) - { - #probs = flavec(ddep,la,mu,K,r,lx,k1) * probs # speciation event - probs = lambdamu(0:(lx - 1) + k1,c(pars1[1:3],r),ddep)[[1]] * probs - } - cp <- check_probs(loglik,probs,verbose); loglik <- cp[[1]]; probs<- cp[[2]]; - } - } else { - probs = rep(0,lx + 1) - probs[1 + missnumspec] = 1 - for(k in (S + 2 - soc):2) - { - k1 = k + (soc - 2) - #y = deSolve::ode(probs,-brts[k:(k-1)],dd_loglik_bw_rhs,c(pars1,k1,ddep),rtol = reltol,atol = abstol,method = methode) - #probs2 = y[2,2:(lx+2)] - probs = dd_loglik_M_bw(pars1,lx,k1,ddep,tt = abs(brts[k] - brts[k-1]),probs[1:lx]) - probs = c(probs,0) - if(k > soc) - { - #probs = c(flavec(ddep,la,mu,K,r,lx,k1-1),1) * probs # speciation event - probs = c(lambdamu(0:(lx - 1) + k1 - 1,pars1,ddep)[[1]],1) * probs - } - cp <- check_probs(loglik,probs[1:lx],verbose); loglik <- cp[[1]]; probs[1:lx] <- cp[[2]]; - } - } - if(probs[1 + (cond != 3) * missnumspec] <= 0 | loglik == -Inf) - { + if((mu == 0 & (ddep == 2 | ddep == 2.1 | ddep == 2.2)) | (la == 0 & (ddep == 4 | ddep == 4.1 | ddep == 4.2)) | (la <= mu)) + { + if(verbose) cat("These parameter values cannot satisfy lambda(N) = mu(N) for a positive and finite N.\n") + loglik = -Inf + } else { + if(((ddep == 1 | ddep == 5) & ceiling(la/(la - mu) * (r + 1) * K) < (S + missnumspec)) | ((ddep == 1.3) & ((S + missnumspec) > ceiling(K)))) + { loglik = -Inf - } else { - loglik = loglik + (cond != 3 | soc == 1) * log(probs[1 + (cond != 3) * missnumspec]) - lgamma(S + missnumspec + 1) + lgamma(S + 1) + lgamma(missnumspec + 1) - - logliknorm = 0 - if(cond == 1 | cond == 2) + } else { + loglik = (btorph == 0) * lgamma(S) + if(cond != 3) { - probsn = rep(0,lx) - probsn[1] = 1 # change if other species at stem or crown age - k = soc - t1 = brts[1] - t2 = brts[S + 2 - soc] - #y = deSolve::ode(probsn,c(t1,t2),rhs_func,c(pars1,k,ddep),rtol = reltol,atol = abstol,method = methode); - #probsn = y[2,2:(lx+1)] - probsn = dd_loglik_M(pars1,lx,k,ddep,tt = abs(t2 - t1),probsn) - if(soc == 1) { aux = 1:lx } - if(soc == 2) { aux = (2:(lx+1)) * (3:(lx+2))/6 } - probsc = probsn/aux - if(cond == 1) { logliknorm = log(sum(probsc)) } - if(cond == 2) { logliknorm = log(probsc[S + missnumspec - soc + 1])} + probs = rep(0,lx) + probs[1] = 1 # change if other species at stem/crown age + for(k in 2:(S + 2 - soc)) + { + k1 = k + (soc - 2) + #y = deSolve::ode(probs,brts[(k-1):k],rhs_func,c(pars1,k1,ddep),rtol = reltolint,atol = abstolint,method = methode) + #probs2 = y[2,2:(lx+1)] + probs = dd_loglik_M(pars1,lx,k1,ddep,tt = abs(brts[k] - brts[k-1]),probs) + if(is.na(sum(probs)) && pars1[2]/pars1[1] < 1E-4 && missnumspec == 0) + { + loglik = dd_loglik_high_lambda(pars1 = pars1,pars2 = pars2,brts = brts) + if(verbose) cat('High lambda approximation has been applied.\n') + return(loglik) + } + if(k < (S + 2 - soc)) + { + #probs = flavec(ddep,la,mu,K,r,lx,k1) * probs # speciation event + probs = lambdamu(0:(lx - 1) + k1,c(pars1[1:3],r),ddep)[[1]] * probs + } + cp <- check_probs(loglik,probs,verbose); loglik <- cp[[1]]; probs<- cp[[2]]; + } + } else { + probs = rep(0,lx + 1) + probs[1 + missnumspec] = 1 + for(k in (S + 2 - soc):2) + { + k1 = k + (soc - 2) + #y = deSolve::ode(probs,-brts[k:(k-1)],dd_loglik_bw_rhs,c(pars1,k1,ddep),rtol = reltolint,atol = abstolint,method = methode) + #probs2 = y[2,2:(lx+2)] + probs = dd_loglik_M_bw(pars1,lx,k1,ddep,tt = abs(brts[k] - brts[k-1]),probs[1:lx]) + probs = c(probs,0) + if(k > soc) + { + #probs = c(flavec(ddep,la,mu,K,r,lx,k1-1),1) * probs # speciation event + probs = c(lambdamu(0:(lx - 1) + k1 - 1,pars1,ddep)[[1]],1) * probs + } + cp <- check_probs(loglik,probs[1:lx],verbose); loglik <- cp[[1]]; probs[1:lx] <- cp[[2]]; + } } - if(cond == 3) - { - #probsn = rep(0,lx + 1) - #probsn[S + missnumspec + 1] = 1 #/ (S + missnumspec) - #TT = max(1,1/abs(la - mu)) * 100000000 * max(abs(brts)) # make this more efficient later - #y = deSolve::ode(probsn,c(0,TT),dd_loglik_bw_rhs,c(pars1,0,ddep),rtol = reltol,atol = abstol,method = methode) - #logliknorm = log(y[2,lx + 2]) - probsn = rep(0,lx + 1) - probsn[2] = 1 - MM = dd_loglik_M_aux(pars1,lx + 1,k = 0,ddep) - MM = MM[-1,-1] - #probsn = SparseM::solve(-MM,probsn[2:(lx + 1)]) - MMinv = SparseM::solve(MM) - probsn = -MMinv %*% probsn[2:(lx + 1)] - logliknorm = log(probsn[S + missnumspec]) - if(soc == 2) - { + if(probs[1 + (cond != 3) * missnumspec] <= 0 | loglik == -Inf) + { + loglik = -Inf + } else { + loglik = loglik + (cond != 3 | soc == 1) * log(probs[1 + (cond != 3) * missnumspec]) - lgamma(S + missnumspec + 1) + lgamma(S + 1) + lgamma(missnumspec + 1) + + logliknorm = 0 + if(cond == 1 | cond == 2) + { + probsn = rep(0,lx) + probsn[1] = 1 # change if other species at stem or crown age + k = soc + t1 = brts[1] + t2 = brts[S + 2 - soc] + #y = deSolve::ode(probsn,c(t1,t2),rhs_func,c(pars1,k,ddep),rtol = reltolint,atol = abstolint,method = methode); + #probsn = y[2,2:(lx+1)] + probsn = dd_loglik_M(pars1,lx,k,ddep,tt = abs(t2 - t1),probsn) + if(soc == 1) { aux = 1:lx } + if(soc == 2) { aux = (2:(lx+1)) * (3:(lx+2))/6 } + probsc = probsn/aux + if(cond == 1) { logliknorm = log(sum(probsc)) } + if(cond == 2) { logliknorm = log(probsc[S + missnumspec - soc + 1])} + } + if(cond == 3) + { + #probsn = rep(0,lx + 1) + #probsn[S + missnumspec + 1] = 1 #/ (S + missnumspec) + #TT = max(1,1/abs(la - mu)) * 100000000 * max(abs(brts)) # make this more efficient later + #y = deSolve::ode(probsn,c(0,TT),dd_loglik_bw_rhs,c(pars1,0,ddep),rtol = reltolint,atol = abstolint,method = methode) + #logliknorm = log(y[2,lx + 2]) + probsn = rep(0,lx + 1) + probsn[2] = 1 + MM = dd_loglik_M_aux(pars1,lx + 1,k = 0,ddep) + MM = MM[-1,-1] + #probsn = SparseM::solve(-MM,probsn[2:(lx + 1)]) + MMinv = SparseM::solve(MM) + probsn = -MMinv %*% probsn[2:(lx + 1)] + logliknorm = log(probsn[S + missnumspec]) + if(soc == 2) + { #probsn = rep(0,lx + 1) #probsn[1:lx] = probs[1:lx] #probsn = c(flavec(ddep,la,mu,K,r,lx,1),1) * probsn # speciation event #probsn = c(lambdamu(0:(lx - 1) + 1,pars1,ddep)[[1]],1) * probsn # speciation event - #y = deSolve::ode(probsn,c(max(abs(brts)),TT),dd_loglik_bw_rhs,c(pars1,1,ddep),rtol = reltol,atol = abstol,method = methode) + #y = deSolve::ode(probsn,c(max(abs(brts)),TT),dd_loglik_bw_rhs,c(pars1,1,ddep),rtol = reltolint,atol = abstolint,method = methode) #logliknorm = logliknorm - log(y[2,lx + 2]) probsn2 = rep(0,lx) probsn2 = lambdamu(0:(lx - 1) + 1,pars1,ddep)[[1]] * probs[1:lx] @@ -466,27 +501,27 @@ if((mu == 0 & (ddep == 2 | ddep == 2.1 | ddep == 2.2)) | (la == 0 & (ddep == 4 | #probsn2 = SparseM::solve(-MM,probsn2[1:lx]) probsn2 = -MMinv %*% probsn2[1:lx] logliknorm = logliknorm - log(probsn2[1]) - } + } + } + loglik = loglik - logliknorm } - loglik = loglik - logliknorm - } + } + }} + if(verbose) + { + s1 = sprintf('Parameters: %f %f %f',pars1[1],pars1[2],pars1[3]) + if(ddep == 5) {s1 = sprintf('%s %f',s1,pars1[4])} + s2 = sprintf(', Loglikelihood: %f',loglik) + cat(s1,s2,"\n",sep = "") + utils::flush.console() } -}} -if(verbose) -{ - s1 = sprintf('Parameters: %f %f %f',pars1[1],pars1[2],pars1[3]) - if(ddep == 5) {s1 = sprintf('%s %f',s1,pars1[4])} - s2 = sprintf(', Loglikelihood: %f',loglik) - cat(s1,s2,"\n",sep = "") - utils::flush.console() -} -} -loglik = as.numeric(loglik) -if(is.nan(loglik) | is.na(loglik) | loglik == Inf) -{ + } + loglik = as.numeric(loglik) + if(is.nan(loglik) | is.na(loglik) | loglik == Inf) + { loglik = -Inf -} -return(loglik) + } + return(loglik) } dd_int <- function(initprobs,tvec,rhs_func,pars,rtol,atol,method) @@ -522,13 +557,13 @@ dd_integrate <- function(initprobs,tvec,rhs_func,pars,rtol,atol,method) #rhs_func_name <- 'no_name' #if(is.character(rhs_func)) #{ - rhs_func_name <- rhs_func - #if(rhs_func_name != 'dd_loglik_rhs' & rhs_func_name != 'dd_loglik_bw_rhs') - #{ - rhs_func = match.fun(rhs_func) - #} + rhs_func_name <- rhs_func + #if(rhs_func_name != 'dd_loglik_rhs' & rhs_func_name != 'dd_loglik_bw_rhs') + #{ + rhs_func = match.fun(rhs_func) + #} #} - if(rhs_func_name == 'dd_loglik_rhs' || rhs_func_name == 'dd_loglik_bw_rhs') + if(rhs_func_name == 'dd_loglik_rhs' || rhs_func_name == 'dd_loglik_log_rhs' || rhs_func_name == 'dd_loglik_bw_rhs') { parsvec = c(dd_loglik_rhs_precomp(pars,initprobs),pars[length(pars) - 1]) } else @@ -549,7 +584,8 @@ dd_integrate <- function(initprobs,tvec,rhs_func,pars,rtol,atol,method) dd_rhs_odeint_map = list( 'dd_loglik_rhs' = dd_integrate_odeint, - 'dd_loglik_bw_rhs' = dd_integrate_bw_odeint + 'dd_loglik_bw_rhs' = dd_integrate_bw_odeint, + 'dd_loglik_log_rhs' = dd_integrate_log_odeint ) diff --git a/R/dd_loglik_M.R b/R/dd_loglik_M.R index d295ce1..d81ed52 100644 --- a/R/dd_loglik_M.R +++ b/R/dd_loglik_M.R @@ -82,7 +82,8 @@ dd_loglik_M = function(pars,lx,k,ddep,tt,p) pars = changepars(pars) MM = dd_loglik_M_aux(pars,lx,k,ddep) #p = expoRkit::expv(x = MM,v = p,t = tt,m = 50L) - p = expm::expAtv(A = MM,v = p,t = tt,m.max = 50)[[1]] + p = try(expm::expAtv(A = MM,v = p,t = tt,m.max = 50)[[1]]) + if(!is.numeric(p)) p <- rep(0,length(p)) return(p) } diff --git a/R/dd_loglik_rhs.R b/R/dd_loglik_rhs.R index b6b080b..8d8b0d1 100644 --- a/R/dd_loglik_rhs.R +++ b/R/dd_loglik_rhs.R @@ -88,4 +88,19 @@ dd_loglik_rhs = function(t,x,parsvec) xx = c(0,x,0) dx = lavec[(2:(lx+1))+kk-1] * nn[(2:(lx+1))+2*kk-1] * xx[(2:(lx+1))-1] + muvec[(2:(lx+1))+kk+1] * nn[(2:(lx+1))+1] * xx[(2:(lx+1))+1] - (lavec[(2:(lx+1))+kk] + muvec[(2:(lx+1))+kk]) * nn[(2:(lx+1))+kk] * xx[2:(lx+1)] return(list(dx)) +} + +dd_loglik_log_rhs = function(t,x,parsvec) +{ + lv = (length(parsvec) - 1)/3 + lavec = parsvec[1:lv] + muvec = parsvec[(lv + 1):(2 * lv)] + nn = parsvec[(2 * lv + 1):(3 * lv)] + kk = parsvec[length(parsvec)] + lx = length(x) + xx = c(-Inf,x,-Inf) + dx = lavec[(2:(lx+1))+kk-1] * nn[(2:(lx+1))+2*kk-1] * exp(xx[(2:(lx+1))-1] - xx[2:(lx+1)]) + + muvec[(2:(lx+1))+kk+1] * nn[(2:(lx+1))+1] * exp(xx[(2:(lx+1))+1] - xx[2:(lx+1)]) - + (lavec[(2:(lx+1))+kk] + muvec[(2:(lx+1))+kk]) * nn[(2:(lx+1))+kk] + return(list(dx)) } \ No newline at end of file diff --git a/R/dd_utils.R b/R/dd_utils.R index 5815c58..3577199 100644 --- a/R/dd_utils.R +++ b/R/dd_utils.R @@ -93,17 +93,17 @@ brts2phylo <- function(times,root=FALSE,tip.label=NULL) #' conv(1:10,1:10) #' #' @export conv -conv = function(x,y) +conv <- function(x,y) { - lx = length(x) - ly = length(y) - lxy = length(x) + length(y) - x = c(x,rep(0,lxy - lx)) - y = c(y,rep(0,lxy - ly)) - cvxy = rep(0,lxy) + lx <- length(x) + ly <- length(y) + lxy <- length(x) + length(y) + x <- c(x,rep(0,lxy - lx)) + y <- c(y,rep(0,lxy - ly)) + cvxy <- rep(0,lxy) for(i in 2:lxy) { - cvxy[i] = crossprod(x[(i-1):1],y[1:(i-1)]) + cvxy[i] <- crossprod(x[(i-1):1],y[1:(i-1)]) } return(cvxy[2:lxy]) } @@ -514,11 +514,13 @@ sample2 = function(x,size,replace = FALSE,prob = NULL) #' #' @param fun Function to be optimized #' @param trparsopt Initial guess of the parameters to be optimized -#' @param ... Any other arguments of the function to be optimimzed, or settings +#' @param ... Any other arguments of the function to be optimimized, or settings #' of the optimization routine -#' @param optimpars Parameters of the optimization: relative tolerance in -#' function arguments, relative tolerance in function value, absolute tolerance -#' in function arguments, and maximum number of iterations +#' @param optimpars Parameters of the optimization: 1) relative tolerance in +#' function arguments, 2) relative tolerance in function value, 3) absolute +#' tolerance in function arguments as well as the function value, 4) +#' maximum number of iterations and 5) TRUE/FALSE flag to allow verbose output, +#' default is TRUE #' @return \item{out}{ A list containing optimal function arguments #' (\code{par}, the optimal function value (\code{fvalues}) and whether the #' optimization converged (\code{conv})}. @@ -529,13 +531,14 @@ sample2 = function(x,size,replace = FALSE,prob = NULL) #' cat("No examples") #' #' @export simplex -simplex = function(fun,trparsopt,optimpars,...) +simplex = function(fun, trparsopt, optimpars, ...) { numpar = length(trparsopt) reltolx = optimpars[1] reltolf = optimpars[2] abstolx = optimpars[3] maxiter = optimpars[4] + if(length(optimpars) > 4) verbose <- optimpars[5] else verbose <- 1 ## Setting up initial simplex v = t(matrix(rep(trparsopt,each = numpar + 1),nrow = numpar + 1)) @@ -557,17 +560,19 @@ simplex = function(fun,trparsopt,optimpars,...) { fv[i] = -fun(trparsopt = v[,i], ...) } - - how = "initial" itercount = 1 - string = itercount - for(i in 1:numpar) - { - string = paste(string, untransform_pars(v[i,1]), sep = " ") + + if(verbose) { + how = "initial" + string = itercount + for(i in 1:numpar) + { + string = paste(string, untransform_pars(v[i,1]), sep = " ") + } + string = paste(string, -fv[1], how, "\n", sep = " ") + cat(string) + #utils::flush.console() } - string = paste(string, -fv[1], how, "\n", sep = " ") - cat(string) - utils::flush.console() tmp = order(fv) if(numpar == 1) @@ -667,21 +672,23 @@ simplex = function(fun,trparsopt,optimpars,...) } fv = fv[tmp] itercount = itercount + 1 - string = itercount; - for(i in 1:numpar) - { + if(verbose) { + string = itercount; + for(i in 1:numpar) + { string = paste(string, untransform_pars(v[i,1]), sep = " ") + } + string = paste(string, -fv[1], how, "\n", sep = " ") + cat(string) } - string = paste(string, -fv[1], how, "\n", sep = " ") - cat(string) - utils::flush.console() + v2 = t(matrix(rep(v[,1],each = numpar + 1),nrow = numpar + 1)) } if(itercount < maxiter) { - cat("Optimization has terminated successfully.","\n") + if (verbose) cat("Optimization has terminated successfully.","\n") } else { - cat("Maximum number of iterations has been exceeded.","\n") + if (verbose) cat("Maximum number of iterations has been exceeded.","\n") } out = list(par = v[,1], fvalues = -fv[1], conv = as.numeric(itercount > maxiter)) invisible(out) @@ -696,18 +703,20 @@ simplex = function(fun,trparsopt,optimpars,...) #' #' @param optimmethod The method to use for optimization, either 'simplex' or #' 'subplex' -#' @param optimpars Parameters of the optimization: relative tolerance in -#' function arguments, relative tolerance in function value, absolute tolerance -#' in function arguments as well as the function value, and maximum number of iterations +#' @param optimpars Parameters of the optimization: 1) relative tolerance in +#' function arguments, 2) relative tolerance in function value, 3) absolute +#' tolerance in function arguments as well as the function value, 4) +#' maximum number of iterations and 5) TRUE/FALSE flag to allow verbose output +#' when using the simplex method, default is TRUE #' @param num_cycles Number of cycles of the optimization. When set to Inf, the #' optimization will be repeated until the result is, within the tolerance, #' equal to the starting values, with a maximum of 10 cycles. #' @param fun Function to be optimized #' @param trparsopt Initial guess of the parameters to be optimized -#' @param jitter Perturbation of an initial parameter value when precisely equal to 0.5; -#' this is only relevant when subplex is chosen. The default value is 0, so no jitter -#' is applied. A recommended value when using it is 1E-5. -#' @param ... Any other arguments of the function to be optimimzed, or settings +#' @param jitter Perturbation of an initial parameter value when precisely equal +#' to 0.5; this is only relevant when subplex is chosen. The default value is +#' 0, so no jitter is applied. A recommended value when using it is 1E-5. +#' @param ... Any other arguments of the function to be optimimized, or settings #' of the optimization routine #' @return \item{out}{ A list containing optimal function arguments #' (\code{par}, the optimal function value (\code{fvalues}) and whether the @@ -782,7 +791,7 @@ optimizer <- function( itermax = optimpars[4], packages = c('DDD')), fun = fun, - + ...))$optim outnew <- list(par = outnew$bestmem, fvalues = -outnew$bestval, conv = 0) } else if(substr(optimmethod,1,7) == 'optim::') @@ -791,14 +800,14 @@ optimizer <- function( { return(-fun(trparsopt = trparsopt, ...)) } - outnew <- suppressWarnings(optim(par = trparsopt, - fn = minfun3, - method = substr(optimmethod,8,nchar(optimmethod)), - control = list(reltol = optimpars[2], - abstol = optimpars[3], - maxit = optimpars[4]), - fun = fun, - ...)) + outnew <- suppressWarnings(stats::optim(par = trparsopt, + fn = minfun3, + method = substr(optimmethod,8,nchar(optimmethod)), + control = list(reltol = optimpars[2], + abstol = optimpars[3], + maxit = optimpars[4]), + fun = fun, + ...)) outnew <- list(par = outnew$par, fvalues = -outnew$value, conv = outnew$convergence) } if(cy > 1 & (any(is.na(outnew$par)) | any(is.nan(outnew$par)) | is.na(outnew$fvalues) | is.nan(outnew$fvalues) | outnew$conv != 0)) @@ -826,6 +835,8 @@ optimizer <- function( } return(out) } + + #' @name transform_pars #' @title Transforming parameters from -Inf to Inf into parameters #' from -1 to 1 diff --git a/README.md b/README.md index b14fdad..b3f1855 100644 --- a/README.md +++ b/README.md @@ -5,7 +5,7 @@ [![](http://cranlogs.r-pkg.org/badges/DDD)](https://CRAN.R-project.org/package=DDD) Branch|[![GitHub Actions logo](man/figures/github_actions_logo.png)](https://github.com/features/actions)|[![Codecov logo](man/figures/Codecov.png)](https://www.codecov.io) ---------|------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|-------------------------------------------------------------------------------------------------------------------------------------------------------------------- +------------------|--------------------------------------------------------------------------------------------------------------------------------------------------------------------------|-------------------------------------------------------------------------------------------------------------------------------------------------------------------- `master`|[![Build Status](https://github.com/rsetienne/DDD/workflows/R-CMD-check/badge.svg?branch=master)](https://github.com/rsetienne/DDD/actions)|[![codecov.io](https://codecov.io/github/rsetienne/DDD/coverage.svg?branch=master)](https://codecov.io/github/rsetienne/DDD/branch/master) `develop`|[![Build Status](https://github.com/rsetienne/DDD/workflows/R-CMD-check/badge.svg?branch=develop)](https://github.com/rsetienne/DDD/actions)|[![codecov.io](https://codecov.io/github/rsetienne/DDD/coverage.svg?branch=develop)](https://codecov.io/github/rsetienne/DDD/branch/develop) diff --git a/man/bd_ML.Rd b/man/bd_ML.Rd index 463e64b..0b800d2 100644 --- a/man/bd_ML.Rd +++ b/man/bd_ML.Rd @@ -20,7 +20,7 @@ bd_ML( tol = c(0.001, 1e-04, 1e-06), maxiter = 1000 * round((1.25)^length(idparsopt)), changeloglikifnoconv = FALSE, - optimmethod = "subplex", + optimmethod = "simplex", num_cycles = 1, methode = "odeint::runge_kutta_cash_karp54", verbose = FALSE @@ -85,7 +85,7 @@ abstolx = absolute tolerance of parameter values in optimization} does not converge} \item{optimmethod}{Method used in optimization of the likelihood. Current -default is 'subplex'. Alternative is 'simplex' (default of previous +default is 'simplex'. Alternative is 'subplex' (default of previous versions)} \item{num_cycles}{the number of cycles of opimization. If set at Inf, it will diff --git a/man/dd_KI_ML.Rd b/man/dd_KI_ML.Rd index 99d9d3e..eb869cb 100644 --- a/man/dd_KI_ML.Rd +++ b/man/dd_KI_ML.Rd @@ -24,7 +24,7 @@ dd_KI_ML( tol = c(0.001, 1e-04, 1e-06), maxiter = 1000 * round((1.25)^length(idparsopt)), changeloglikifnoconv = FALSE, - optimmethod = "subplex", + optimmethod = "simplex", num_cycles = 1, methode = "analytical", correction = TRUE, @@ -110,7 +110,7 @@ abstolx = absolute tolerance of parameter values in optimization} does not converge} \item{optimmethod}{Method used in optimization of the likelihood. Current -default is 'subplex'. Alternative is 'simplex' (default of previous +default is 'simplex'. Alternative is 'subplex' (default of previous versions)} \item{num_cycles}{the number of cycles of opimization. If set at Inf, it will diff --git a/man/dd_LR.Rd b/man/dd_LR.Rd index f45996f..1f3b76d 100644 --- a/man/dd_LR.Rd +++ b/man/dd_LR.Rd @@ -22,7 +22,7 @@ dd_LR( tol = c(0.001, 1e-04, 1e-06), maxiter = 2000, changeloglikifnoconv = FALSE, - optimmethod = "subplex", + optimmethod = "simplex", methode = "analytical" ) } @@ -96,7 +96,7 @@ tolerance of parameter values in optimization} does not converge} \item{optimmethod}{Method used in optimization of the likelihood. Current -default is 'subplex'. Alternative is 'simplex' (default of previous +default is 'simplex'. Alternative is 'subplex' (default of previous versions)} \item{methode}{The method used to solve the master equation, default is diff --git a/man/dd_ML.Rd b/man/dd_ML.Rd index 5505077..198514a 100644 --- a/man/dd_ML.Rd +++ b/man/dd_ML.Rd @@ -18,9 +18,11 @@ dd_ML( btorph = 1, soc = 2, tol = c(0.001, 1e-04, 1e-06), + tolint = c(1e-10, 1e-08), + probs_threshold = 0, maxiter = 1000 * round((1.25)^length(idparsopt)), changeloglikifnoconv = FALSE, - optimmethod = "subplex", + optimmethod = "simplex", num_cycles = 1, methode = "analytical", verbose = FALSE @@ -95,10 +97,17 @@ the phylogeny (1)} \item{soc}{Sets whether stem or crown age should be used (1 or 2)} -\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{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{tolint}{Sets the tolerance of the numerical integration. Consists of: \cr +absoltint = absolute tolerance and \cr +reltolint = relative tolerance.} + +\item{probs_threshold}{Sets the threshold of the probability below which +logarithmic integration must be used. Default is 0.} \item{maxiter}{Sets the maximum number of iterations in the optimization} @@ -106,7 +115,7 @@ tolerance of parameter values in optimization} does not converge} \item{optimmethod}{Method used in optimization of the likelihood. Current -default is 'subplex'. Alternative is 'simplex' (default of previous +default is 'simplex'. Alternative is 'subplex' (default of previous versions)} \item{num_cycles}{the number of cycles of opimization. If set at Inf, it will diff --git a/man/dd_MS_ML.Rd b/man/dd_MS_ML.Rd index 97e95c0..e66fe7e 100644 --- a/man/dd_MS_ML.Rd +++ b/man/dd_MS_ML.Rd @@ -24,7 +24,7 @@ dd_MS_ML( tol = c(0.001, 1e-04, 1e-06), maxiter = 1000 * round((1.25)^length(idparsopt)), changeloglikifnoconv = FALSE, - optimmethod = "subplex", + optimmethod = "simplex", num_cycles = 1, methode = "ode45", correction = FALSE, @@ -108,7 +108,7 @@ tolerance of parameter values in optimization} does not converge} \item{optimmethod}{Method used in optimization of the likelihood. Current -default is 'subplex'. Alternative is 'simplex' (default of previous +default is 'simplex'. Alternative is 'subplex' (default of previous versions)} \item{num_cycles}{the number of cycles of opimization. If set at Inf, it will diff --git a/man/dd_SR_ML.Rd b/man/dd_SR_ML.Rd index d009c40..1c8d06f 100644 --- a/man/dd_SR_ML.Rd +++ b/man/dd_SR_ML.Rd @@ -23,7 +23,7 @@ dd_SR_ML( tol = c(0.001, 1e-04, 1e-06), maxiter = 1000 * round((1.25)^length(idparsopt)), changeloglikifnoconv = FALSE, - optimmethod = "subplex", + optimmethod = "simplex", num_cycles = 1, methode = "analytical", verbose = FALSE @@ -105,7 +105,7 @@ tolerance of parameter values in optimization} does not converge} \item{optimmethod}{Method used in optimization of the likelihood. Current -default is 'subplex'. Alternative is 'simplex' (default of previous +default is 'simplex'. Alternative is 'subplex' (default of previous versions)} \item{num_cycles}{the number of cycles of opimization. If set at Inf, it will diff --git a/man/dd_loglik.Rd b/man/dd_loglik.Rd index ff85205..bff7b95 100644 --- a/man/dd_loglik.Rd +++ b/man/dd_loglik.Rd @@ -65,7 +65,12 @@ the likelihood is for the branching times (0) or the phylogeny (1) \cr \cr \code{pars2[5]} sets whether the parameters and likelihood should be shown on screen (1) or not (0) \cr \cr \code{pars2[6]} sets whether the first data -point is stem age (1) or crown age (2)} +point is stem age (1) or crown age (2) +\cr \cr \code{pars2[7]} sets the absolute tolerance of the integration +\cr \cr \code{pars2[8]} sets the relative tolerance of the integration +\cr \cr \code{pars2[9]} sets threshold on the number of species for which +logarithmic integration will be used. If the last three elements are not set +ddefault values will be used of 1E-10, 1E-8 and Inf.} \item{brts}{A set of branching times of a phylogeny, all positive} diff --git a/man/optimizer.Rd b/man/optimizer.Rd index d3c2a42..0ef3346 100644 --- a/man/optimizer.Rd +++ b/man/optimizer.Rd @@ -18,9 +18,11 @@ optimizer( \item{optimmethod}{The method to use for optimization, either 'simplex' or 'subplex'} -\item{optimpars}{Parameters of the optimization: relative tolerance in -function arguments, relative tolerance in function value, absolute tolerance -in function arguments as well as the function value, and maximum number of iterations} +\item{optimpars}{Parameters of the optimization: 1) relative tolerance in +function arguments, 2) relative tolerance in function value, 3) absolute +tolerance in function arguments as well as the function value, 4) +maximum number of iterations and 5) TRUE/FALSE flag to allow verbose output +when using the simplex method, default is TRUE} \item{num_cycles}{Number of cycles of the optimization. When set to Inf, the optimization will be repeated until the result is, within the tolerance, @@ -30,11 +32,11 @@ equal to the starting values, with a maximum of 10 cycles.} \item{trparsopt}{Initial guess of the parameters to be optimized} -\item{jitter}{Perturbation of an initial parameter value when precisely equal to 0.5; -this is only relevant when subplex is chosen. The default value is 0, so no jitter -is applied. A recommended value when using it is 1E-5.} +\item{jitter}{Perturbation of an initial parameter value when precisely equal +to 0.5; this is only relevant when subplex is chosen. The default value is +0, so no jitter is applied. A recommended value when using it is 1E-5.} -\item{...}{Any other arguments of the function to be optimimzed, or settings +\item{...}{Any other arguments of the function to be optimimized, or settings of the optimization routine} } \value{ diff --git a/man/simplex.Rd b/man/simplex.Rd index 53842a6..ff0f6ef 100644 --- a/man/simplex.Rd +++ b/man/simplex.Rd @@ -11,11 +11,13 @@ simplex(fun, trparsopt, optimpars, ...) \item{trparsopt}{Initial guess of the parameters to be optimized} -\item{optimpars}{Parameters of the optimization: relative tolerance in -function arguments, relative tolerance in function value, absolute tolerance -in function arguments, and maximum number of iterations} +\item{optimpars}{Parameters of the optimization: 1) relative tolerance in +function arguments, 2) relative tolerance in function value, 3) absolute +tolerance in function arguments as well as the function value, 4) +maximum number of iterations and 5) TRUE/FALSE flag to allow verbose output, +default is TRUE} -\item{...}{Any other arguments of the function to be optimimzed, or settings +\item{...}{Any other arguments of the function to be optimimized, or settings of the optimization routine} } \value{ diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 5376c6d..17f20ff 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -27,6 +27,22 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } +// dd_integrate_log_odeint +NumericVector dd_integrate_log_odeint(NumericVector ry, NumericVector times, NumericVector pars, double atol, double rtol, std::string stepper); +RcppExport SEXP _DDD_dd_integrate_log_odeint(SEXP rySEXP, SEXP timesSEXP, SEXP parsSEXP, SEXP atolSEXP, SEXP rtolSEXP, SEXP stepperSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< NumericVector >::type ry(rySEXP); + Rcpp::traits::input_parameter< NumericVector >::type times(timesSEXP); + Rcpp::traits::input_parameter< NumericVector >::type pars(parsSEXP); + Rcpp::traits::input_parameter< double >::type atol(atolSEXP); + Rcpp::traits::input_parameter< double >::type rtol(rtolSEXP); + Rcpp::traits::input_parameter< std::string >::type stepper(stepperSEXP); + rcpp_result_gen = Rcpp::wrap(dd_integrate_log_odeint(ry, times, pars, atol, rtol, stepper)); + return rcpp_result_gen; +END_RCPP +} // dd_integrate_odeint NumericVector dd_integrate_odeint(NumericVector ry, NumericVector times, NumericVector pars, double atol, double rtol, std::string stepper); RcppExport SEXP _DDD_dd_integrate_odeint(SEXP rySEXP, SEXP timesSEXP, SEXP parsSEXP, SEXP atolSEXP, SEXP rtolSEXP, SEXP stepperSEXP) { @@ -94,6 +110,7 @@ END_RCPP static const R_CallMethodDef CallEntries[] = { {"_DDD_dd_integrate_bw_odeint", (DL_FUNC) &_DDD_dd_integrate_bw_odeint, 6}, + {"_DDD_dd_integrate_log_odeint", (DL_FUNC) &_DDD_dd_integrate_log_odeint, 6}, {"_DDD_dd_integrate_odeint", (DL_FUNC) &_DDD_dd_integrate_odeint, 6}, {"_DDD_dd_integrate_td_odeint", (DL_FUNC) &_DDD_dd_integrate_td_odeint, 6}, {"_DDD_dd_logliknorm1_odeint", (DL_FUNC) &_DDD_dd_logliknorm1_odeint, 6}, diff --git a/src/dd_integrate_log_odeint.cpp b/src/dd_integrate_log_odeint.cpp new file mode 100644 index 0000000..984eec3 --- /dev/null +++ b/src/dd_integrate_log_odeint.cpp @@ -0,0 +1,79 @@ +//' @useDynLib DDD + + +#define STRICT_R_HEADERS +#include "config.h" +#include +#include +#include +#include "odeint_helper.h" +#include + +using namespace Rcpp; + + +class ode_log_rhs +{ +public: + ode_log_rhs(NumericVector parsvec) + { + const size_t lv = (parsvec.size() - 1) / 3; + lavec.resize(lv, 0); + muvec.resize(lv, 0); + nn.resize(lv, 0); + for (size_t i = 0; i < lv; ++i) { + lavec[i] = parsvec[i]; // parsvec[1:lv] + muvec[i] = parsvec[lv + i]; // parsvec[(lv + 1):(2 * lv)] + nn[i] = parsvec[2 * lv + i]; // parsvec[(2 * lv + 1):(3 * lv)] + } + kk = static_cast(parsvec[parsvec.size() - 1]); + } + + void operator()(const std::vector& xx, std::vector& dx, double /* t */) + { + // R code: + // lx = length(x) + // xx = c(0,x,0) + // dx = lavec[(2:(lx+1))+kk-1] * nn[(2:(lx+1))+2*kk-1] * exp(xx[(2:(lx+1))-1] - xx[2:(lx+1)]) + // + muvec[(2:(lx+1))+kk+1] * nn[(2:(lx+1))+1] * exp(xx[(2:(lx+1))+1] - xx[2:(lx+1)]) + // - (lavec[(2:(lx+1))+kk] + muvec[(2:(lx+1))+kk]) * nn[(2:(lx+1))+kk] + // return list(dx) + + dx.front() = dx.back() = 0.0; + const size_t lx = xx.size() - 1; + dx[1] = muvec[2 + kk] * nn[2] * exp(xx[2] - xx[1]) + - (lavec[1 + kk] + muvec[1 + kk]) * nn[1 + kk]; + dx[lx - 1] = lavec[lx - 2 + kk] * nn[lx - 2 + 2*kk] * exp(xx[lx - 2] - xx[lx - 1]) + - (lavec[lx - 1 + kk] + muvec[lx - 1 + kk]) * nn[lx - 1 + kk]; + for (size_t i = 2; i < lx - 1; ++i) { + const size_t i0 = i - 1; + const size_t i1 = i + 1; + dx[i] = lavec[i0 + kk] * nn[i0 + 2*kk] * exp(xx[i0] - xx[i]) + + muvec[i1 + kk] * nn[i1] * exp(xx[i1] - xx[i]) + - (lavec[i + kk] + muvec[i + kk]) * nn[i + kk]; + } + } + +private: + size_t kk; + std::vector lavec; + std::vector muvec; + std::vector nn; +}; + + +// [[Rcpp::export]] +NumericVector dd_integrate_log_odeint(NumericVector ry, + NumericVector times, + NumericVector pars, + double atol, + double rtol, + std::string stepper) +{ + std::vector y(ry.size() + 2, 0.0); // [0,y,0] + std::copy(ry.begin(), ry.end(), y.begin() + 1); + + auto rhs_obj = ode_log_rhs(pars); + odeint_helper::integrate(stepper, std::ref(rhs_obj), y, times[0], times[1], 0.1 * (times[1] - times[0]), atol, rtol); + return NumericVector(y.cbegin() + 1, y.cend() - 1); +}