diff --git a/DDD.Rproj b/DDD.Rproj index b762949..301db4c 100644 --- a/DDD.Rproj +++ b/DDD.Rproj @@ -1,5 +1,4 @@ Version: 1.0 -ProjectId: cbff7951-4a06-4670-bbdf-df0ae63c05bf RestoreWorkspace: Default SaveWorkspace: Default 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/dd_ML.R b/R/dd_ML.R index 44541d7..e9cc9df 100644 --- a/R/dd_ML.R +++ b/R/dd_ML.R @@ -97,6 +97,8 @@ parsfixdefault = function(ddmodel,brts,missnumspec,idparsopt) #' 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 k_threshold Sets the threshold of number of species above which +#' logarithmic integration must be used. Default is Inf. #' @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 @@ -148,6 +150,7 @@ dd_ML = function( soc = 2, tol = c(1E-3, 1E-4, 1E-6), tolint = c(1E-10,1E-8), + k_threshold = 3, maxiter = 1000 * round((1.25)^length(idparsopt)), changeloglikifnoconv = FALSE, optimmethod = 'simplex', @@ -186,7 +189,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,abstolint = tolint[1],reltolint = tolint[2]) + pars2 = c(res,ddmodel,cond,btorph,verbose,soc,tol,maxiter,abstolint = tolint[1],reltolint = tolint[2],k_threshold = k_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_loglik.R b/R/dd_loglik.R index 08ee548..c1bf988 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') @@ -183,6 +188,7 @@ dd_loglik1 = function(pars1,pars2,brts,missnumspec,methode = 'odeint::runge_kutt } else { 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['k_threshold'])) k_threshold <- Inf else k_threshold <- pars2['k_threshold'] brts = -sort(abs(as.numeric(brts)),decreasing = TRUE) if(sum(brts == 0) == 0) { @@ -209,9 +215,20 @@ 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 for(k in 2:(S + 2 - soc)) { k1 = k + (soc - 2) + probs <- probs[1:lx] + if(k >= k_threshold) { + lx <- min(lx,which(probs == 0) - 1) + probs <- probs[1:lx] + if(k == k_threshold) { + 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(k >= k_threshold) { + 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) @@ -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((k < k_threshold & 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(k < k_threshold) { + 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) @@ -315,146 +348,146 @@ 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 { -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 -} 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) { - loglik = -Inf + brts[length(brts) + 1] = 0 + } + S = length(brts) + (soc - 2) + if(min(pars1) < 0) + { + 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 = 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(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 = 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])} + 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 = 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) - { + 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 @@ -468,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) @@ -524,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 @@ -551,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/man/dd_ML.Rd b/man/dd_ML.Rd index eaa79cc..4da3b84 100644 --- a/man/dd_ML.Rd +++ b/man/dd_ML.Rd @@ -19,6 +19,7 @@ dd_ML( soc = 2, tol = c(0.001, 1e-04, 1e-06), tolint = c(1e-10, 1e-08), + k_threshold = 3, maxiter = 1000 * round((1.25)^length(idparsopt)), changeloglikifnoconv = FALSE, optimmethod = "simplex", @@ -104,6 +105,9 @@ 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{k_threshold}{Sets the threshold of number of species above which +logarithmic integration must be used. Default is Inf.} + \item{maxiter}{Sets the maximum number of iterations in the optimization} \item{changeloglikifnoconv}{if TRUE the loglik will be set to -Inf if ML 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/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); +}