From d7e3dc62edd910f9717b7f81de61e7150803f1bf Mon Sep 17 00:00:00 2001 From: Thijs Janzen Date: Tue, 4 Mar 2025 15:14:40 +0100 Subject: [PATCH 01/29] add verbose flag to simplex optimization --- DDD.Rproj | 1 + DESCRIPTION | 2 +- R/dd_utils.R | 25 +++++++++++++++---------- man/optimizer.Rd | 3 +++ man/simplex.Rd | 2 +- 5 files changed, 21 insertions(+), 12 deletions(-) diff --git a/DDD.Rproj b/DDD.Rproj index 301db4c..b762949 100644 --- a/DDD.Rproj +++ b/DDD.Rproj @@ -1,4 +1,5 @@ Version: 1.0 +ProjectId: cbff7951-4a06-4670-bbdf-df0ae63c05bf RestoreWorkspace: Default SaveWorkspace: Default diff --git a/DESCRIPTION b/DESCRIPTION index 0a22769..3d5e496 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: DDD Type: Package Title: Diversity-Dependent Diversification -Version: 5.2.3 +Version: 5.2.4 Date: 2024-11-26 Depends: R (>= 3.5.0) Imports: diff --git a/R/dd_utils.R b/R/dd_utils.R index 5815c58..514ea6a 100644 --- a/R/dd_utils.R +++ b/R/dd_utils.R @@ -529,7 +529,7 @@ sample2 = function(x,size,replace = FALSE,prob = NULL) #' cat("No examples") #' #' @export simplex -simplex = function(fun,trparsopt,optimpars,...) +simplex = function(fun,trparsopt,optimpars,verbose = TRUE,...) { numpar = length(trparsopt) reltolx = optimpars[1] @@ -566,7 +566,7 @@ simplex = function(fun,trparsopt,optimpars,...) string = paste(string, untransform_pars(v[i,1]), sep = " ") } string = paste(string, -fv[1], how, "\n", sep = " ") - cat(string) + if (verbose) cat(string) utils::flush.console() tmp = order(fv) @@ -673,15 +673,15 @@ simplex = function(fun,trparsopt,optimpars,...) string = paste(string, untransform_pars(v[i,1]), sep = " ") } string = paste(string, -fv[1], how, "\n", sep = " ") - cat(string) + if (verbose) 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) @@ -707,6 +707,7 @@ simplex = function(fun,trparsopt,optimpars,...) #' @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 verbose if TRUE, prints intermediate output #' @param ... Any other arguments of the function to be optimimzed, or settings #' of the optimization routine #' @return \item{out}{ A list containing optimal function arguments @@ -726,6 +727,7 @@ optimizer <- function( fun, trparsopt, jitter = 0, + verbose = TRUE, ...) { if(num_cycles == Inf) @@ -744,16 +746,19 @@ optimizer <- function( out <- NULL while(cy <= max_cycles) { - if(max_cycles > 1) cat(paste('Cycle ',cy,'\n',sep ='')) + if(max_cycles > 1) { + if (verbose) cat(paste('Cycle ',cy,'\n',sep ='')) + } if(optimmethod == 'simplex') { outnew <- suppressWarnings(simplex(fun = fun, trparsopt = trparsopt, optimpars = optimpars, + verbose = verbose, ...)) } else if(optimmethod == 'subplex') { - minfun1 <- function(fun,trparsopt,...) + minfun1 <- function(fun, trparsopt,...) { return(-fun(trparsopt = trparsopt,...)) } @@ -803,7 +808,7 @@ optimizer <- function( } if(cy > 1 & (any(is.na(outnew$par)) | any(is.nan(outnew$par)) | is.na(outnew$fvalues) | is.nan(outnew$fvalues) | outnew$conv != 0)) { - cat('The last cycle failed; second last cycle result is returned.\n') + if (verbose) cat('The last cycle failed; second last cycle result is returned.\n') return(out) } else { @@ -815,11 +820,11 @@ optimizer <- function( { if(abs(fvalue[cy] - fvalue[cy - 1]) < optimpars[3]) { - if(cy < max_cycles) cat('No more cycles needed.\n') + if(cy < max_cycles) if (verbose) cat('No more cycles needed.\n') cy <- max_cycles } else if(cy == max_cycles) { - cat('More cycles in optimization recommended.\n') + if (verbose) cat('More cycles in optimization recommended.\n') } } cy <- cy + 1 diff --git a/man/optimizer.Rd b/man/optimizer.Rd index d3c2a42..7791aa7 100644 --- a/man/optimizer.Rd +++ b/man/optimizer.Rd @@ -11,6 +11,7 @@ optimizer( fun, trparsopt, jitter = 0, + verbose = TRUE, ... ) } @@ -34,6 +35,8 @@ equal to the starting values, with a maximum of 10 cycles.} 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{verbose}{if TRUE, prints intermediate output} + \item{...}{Any other arguments of the function to be optimimzed, or settings of the optimization routine} } diff --git a/man/simplex.Rd b/man/simplex.Rd index 53842a6..b7d95cf 100644 --- a/man/simplex.Rd +++ b/man/simplex.Rd @@ -4,7 +4,7 @@ \alias{simplex} \title{Carries out optimization using a simplex algorithm (finding a minimum)} \usage{ -simplex(fun, trparsopt, optimpars, ...) +simplex(fun, trparsopt, optimpars, verbose = TRUE, ...) } \arguments{ \item{fun}{Function to be optimized} From bf0e747e744a3ae2317d5ecb6c530a66189bd44c Mon Sep 17 00:00:00 2001 From: Thijs Janzen Date: Tue, 4 Mar 2025 16:36:04 +0100 Subject: [PATCH 02/29] update --- R/bd_ML.R | 4 ++-- R/dd_KI_ML.R | 4 ++-- R/dd_LR.R | 4 ++-- R/dd_ML.R | 4 ++-- R/dd_MS_ML.R | 4 ++-- R/dd_SR_ML.R | 4 ++-- R/dd_utils.R | 6 ++++-- man/bd_ML.Rd | 4 ++-- man/dd_KI_ML.Rd | 4 ++-- man/dd_LR.Rd | 4 ++-- man/dd_ML.Rd | 4 ++-- man/dd_MS_ML.Rd | 4 ++-- man/dd_SR_ML.Rd | 4 ++-- man/simplex.Rd | 4 +++- 14 files changed, 31 insertions(+), 27 deletions(-) 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..e11240b 100644 --- a/R/dd_ML.R +++ b/R/dd_ML.R @@ -99,7 +99,7 @@ parsfixdefault = function(ddmodel,brts,missnumspec,idparsopt) #' @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. @@ -147,7 +147,7 @@ dd_ML = function( 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_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_utils.R b/R/dd_utils.R index 514ea6a..6439a24 100644 --- a/R/dd_utils.R +++ b/R/dd_utils.R @@ -514,11 +514,12 @@ 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 verbose adds verbose output of intermediate evaluations #' @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,7 +530,8 @@ sample2 = function(x,size,replace = FALSE,prob = NULL) #' cat("No examples") #' #' @export simplex -simplex = function(fun,trparsopt,optimpars,verbose = TRUE,...) +simplex = function(fun,trparsopt,optimpars, + verbose = TRUE,...) { numpar = length(trparsopt) reltolx = optimpars[1] 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..801eb85 100644 --- a/man/dd_ML.Rd +++ b/man/dd_ML.Rd @@ -20,7 +20,7 @@ dd_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 @@ -106,7 +106,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/simplex.Rd b/man/simplex.Rd index b7d95cf..15350c4 100644 --- a/man/simplex.Rd +++ b/man/simplex.Rd @@ -15,7 +15,9 @@ simplex(fun, trparsopt, optimpars, verbose = TRUE, ...) function arguments, relative tolerance in function value, absolute tolerance in function arguments, and maximum number of iterations} -\item{...}{Any other arguments of the function to be optimimzed, or settings +\item{verbose}{adds verbose output of intermediate evaluations} + +\item{...}{Any other arguments of the function to be optimimized, or settings of the optimization routine} } \value{ From 7f8e2dad54b12a82c98e713f0a28d544c2f05ae2 Mon Sep 17 00:00:00 2001 From: rsetienne Date: Wed, 5 Mar 2025 13:31:16 +0100 Subject: [PATCH 03/29] update README --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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) From ed2e451d1b2cf1623506177b19e502607b5d7f13 Mon Sep 17 00:00:00 2001 From: rsetienne Date: Wed, 5 Mar 2025 13:38:11 +0100 Subject: [PATCH 04/29] Add verbose option --- R/dd_utils.R | 37 ++++++++++++++++++++----------------- 1 file changed, 20 insertions(+), 17 deletions(-) diff --git a/R/dd_utils.R b/R/dd_utils.R index 5815c58..158195f 100644 --- a/R/dd_utils.R +++ b/R/dd_utils.R @@ -536,6 +536,7 @@ simplex = function(fun,trparsopt,optimpars,...) reltolf = optimpars[2] abstolx = optimpars[3] maxiter = optimpars[4] + verbose <- optimpars[5] ## Setting up initial simplex v = t(matrix(rep(trparsopt,each = numpar + 1),nrow = numpar + 1)) @@ -557,18 +558,18 @@ 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" + itercount = 1 + 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,14 +668,16 @@ 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) + utils::flush.console() } - 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) From 73b3a8edf1762f7ff801d43f8b5fbad81345c236 Mon Sep 17 00:00:00 2001 From: rsetienne Date: Wed, 5 Mar 2025 13:39:35 +0100 Subject: [PATCH 05/29] Version number update --- DESCRIPTION | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 0a22769..fb1fccd 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,8 +1,8 @@ Package: DDD Type: Package Title: Diversity-Dependent Diversification -Version: 5.2.3 -Date: 2024-11-26 +Version: 5.2.4 +Date: 2024-03-04 Depends: R (>= 3.5.0) Imports: deSolve, From 95c84809cf870e4341df8402a6cb9f887fb995c5 Mon Sep 17 00:00:00 2001 From: Thijs Janzen Date: Wed, 5 Mar 2025 14:26:52 +0100 Subject: [PATCH 06/29] make it a verbose flag again --- R/dd_utils.R | 8 ++++---- man/optimizer.Rd | 4 ++-- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/R/dd_utils.R b/R/dd_utils.R index 6439a24..759f571 100644 --- a/R/dd_utils.R +++ b/R/dd_utils.R @@ -709,7 +709,7 @@ simplex = function(fun,trparsopt,optimpars, #' @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 verbose if TRUE, prints intermediate output +#' @param verbose if TRUE, prints intermediate output when using simplex #' @param ... Any other arguments of the function to be optimimzed, or settings #' of the optimization routine #' @return \item{out}{ A list containing optimal function arguments @@ -729,7 +729,7 @@ optimizer <- function( fun, trparsopt, jitter = 0, - verbose = TRUE, + verbose = FALSE, ...) { if(num_cycles == Inf) @@ -810,7 +810,7 @@ optimizer <- function( } if(cy > 1 & (any(is.na(outnew$par)) | any(is.nan(outnew$par)) | is.na(outnew$fvalues) | is.nan(outnew$fvalues) | outnew$conv != 0)) { - if (verbose) cat('The last cycle failed; second last cycle result is returned.\n') + cat('The last cycle failed; second last cycle result is returned.\n') return(out) } else { @@ -826,7 +826,7 @@ optimizer <- function( cy <- max_cycles } else if(cy == max_cycles) { - if (verbose) cat('More cycles in optimization recommended.\n') + cat('More cycles in optimization recommended.\n') } } cy <- cy + 1 diff --git a/man/optimizer.Rd b/man/optimizer.Rd index 7791aa7..58553c0 100644 --- a/man/optimizer.Rd +++ b/man/optimizer.Rd @@ -11,7 +11,7 @@ optimizer( fun, trparsopt, jitter = 0, - verbose = TRUE, + verbose = FALSE, ... ) } @@ -35,7 +35,7 @@ equal to the starting values, with a maximum of 10 cycles.} 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{verbose}{if TRUE, prints intermediate output} +\item{verbose}{if TRUE, prints intermediate output when using simplex} \item{...}{Any other arguments of the function to be optimimzed, or settings of the optimization routine} From 75ea11d5e3e4adc36c24f5534be6f3ecadf2ba12 Mon Sep 17 00:00:00 2001 From: rsetienne Date: Thu, 6 Mar 2025 16:32:47 +0100 Subject: [PATCH 07/29] Add fifth element of optimpars --- R/dd_utils.R | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/R/dd_utils.R b/R/dd_utils.R index 158195f..0610c24 100644 --- a/R/dd_utils.R +++ b/R/dd_utils.R @@ -518,7 +518,7 @@ sample2 = function(x,size,replace = FALSE,prob = NULL) #' 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 +#' in function arguments, maximum number of iterations, and the level of verbosity #' @return \item{out}{ A list containing optimal function arguments #' (\code{par}, the optimal function value (\code{fvalues}) and whether the #' optimization converged (\code{conv})}. @@ -536,7 +536,7 @@ simplex = function(fun,trparsopt,optimpars,...) reltolf = optimpars[2] abstolx = optimpars[3] maxiter = optimpars[4] - verbose <- optimpars[5] + if(length(optimpars) > 4) verbose <- optimpars[5] ## Setting up initial simplex v = t(matrix(rep(trparsopt,each = numpar + 1),nrow = numpar + 1)) From 48b1a30602c983b9da8ed03c46b4597efadd7f08 Mon Sep 17 00:00:00 2001 From: rsetienne Date: Fri, 7 Mar 2025 09:14:32 +0100 Subject: [PATCH 08/29] verbose simplex once more --- R/dd_utils.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/dd_utils.R b/R/dd_utils.R index 0610c24..93429b8 100644 --- a/R/dd_utils.R +++ b/R/dd_utils.R @@ -536,7 +536,7 @@ simplex = function(fun,trparsopt,optimpars,...) reltolf = optimpars[2] abstolx = optimpars[3] maxiter = optimpars[4] - if(length(optimpars) > 4) verbose <- optimpars[5] + 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)) From c3e6389fd406ea0d92fbc8e523904a14d51fd28a Mon Sep 17 00:00:00 2001 From: Thijs Janzen Date: Thu, 13 Mar 2025 14:57:49 +0100 Subject: [PATCH 09/29] Update test-coverage.yaml --- .github/workflows/test-coverage.yaml | 39 ++++++++++++++++++---------- 1 file changed, 26 insertions(+), 13 deletions(-) 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 From c92ac3963b06804775132ae0fbb10bdec9af54db Mon Sep 17 00:00:00 2001 From: Thijs Janzen Date: Thu, 27 Mar 2025 15:35:14 +0100 Subject: [PATCH 10/29] fix simplex bug when verbose = TRUE --- R/dd_utils.R | 56 +++++++++++++++++++++++++++------------------------- 1 file changed, 29 insertions(+), 27 deletions(-) diff --git a/R/dd_utils.R b/R/dd_utils.R index 8ea4b6c..f1b0559 100644 --- a/R/dd_utils.R +++ b/R/dd_utils.R @@ -519,8 +519,6 @@ sample2 = function(x,size,replace = FALSE,prob = NULL) #' @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 verbose adds verbose output of intermediate evaluations -#' in function arguments, maximum number of iterations, and the level of verbosity #' @return \item{out}{ A list containing optimal function arguments #' (\code{par}, the optimal function value (\code{fvalues}) and whether the #' optimization converged (\code{conv})}. @@ -531,8 +529,7 @@ sample2 = function(x,size,replace = FALSE,prob = NULL) #' cat("No examples") #' #' @export simplex -simplex = function(fun,trparsopt,optimpars, - verbose = TRUE,...) +simplex = function(fun, trparsopt, optimpars, ...) { numpar = length(trparsopt) reltolx = optimpars[1] @@ -561,9 +558,10 @@ simplex = function(fun,trparsopt,optimpars, { fv[i] = -fun(trparsopt = v[,i], ...) } + itercount = 1 + if(verbose) { how = "initial" - itercount = 1 string = itercount for(i in 1:numpar) { @@ -574,9 +572,12 @@ simplex = function(fun,trparsopt,optimpars, utils::flush.console() } - string = paste(string, -fv[1], how, "\n", sep = " ") - if (verbose) cat(string) - utils::flush.console() + + if (verbose) { + string = paste(string, -fv[1], how, "\n", sep = " ") + cat(string) + utils::flush.console() + } tmp = order(fv) if(numpar == 1) @@ -687,9 +688,12 @@ simplex = function(fun,trparsopt,optimpars, utils::flush.console() } - string = paste(string, -fv[1], how, "\n", sep = " ") - if (verbose) cat(string) - utils::flush.console() + + if (verbose) { + 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)) } @@ -712,19 +716,19 @@ 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 #' @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 verbose if TRUE, prints intermediate output when using simplex -#' @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 @@ -743,7 +747,6 @@ optimizer <- function( fun, trparsopt, jitter = 0, - verbose = FALSE, ...) { if(num_cycles == Inf) @@ -762,19 +765,16 @@ optimizer <- function( out <- NULL while(cy <= max_cycles) { - if(max_cycles > 1) { - if (verbose) cat(paste('Cycle ',cy,'\n',sep ='')) - } + if(max_cycles > 1) cat(paste('Cycle ',cy,'\n',sep ='')) if(optimmethod == 'simplex') { outnew <- suppressWarnings(simplex(fun = fun, trparsopt = trparsopt, optimpars = optimpars, - verbose = verbose, ...)) } else if(optimmethod == 'subplex') { - minfun1 <- function(fun, trparsopt,...) + minfun1 <- function(fun,trparsopt,...) { return(-fun(trparsopt = trparsopt,...)) } @@ -803,7 +803,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::') @@ -836,7 +836,7 @@ optimizer <- function( { if(abs(fvalue[cy] - fvalue[cy - 1]) < optimpars[3]) { - if(cy < max_cycles) if (verbose) cat('No more cycles needed.\n') + if(cy < max_cycles) cat('No more cycles needed.\n') cy <- max_cycles } else if(cy == max_cycles) { @@ -847,6 +847,8 @@ optimizer <- function( } return(out) } + + #' @name transform_pars #' @title Transforming parameters from -Inf to Inf into parameters #' from -1 to 1 From 0daed3ffd43291cc3319f2b150d43342d7cbbded Mon Sep 17 00:00:00 2001 From: Thijs Janzen Date: Thu, 27 Mar 2025 15:43:25 +0100 Subject: [PATCH 11/29] update documentation --- R/dd_utils.R | 9 ++++++--- man/optimizer.Rd | 19 +++++++++---------- man/simplex.Rd | 12 ++++++------ 3 files changed, 21 insertions(+), 19 deletions(-) diff --git a/R/dd_utils.R b/R/dd_utils.R index f1b0559..a3e7b76 100644 --- a/R/dd_utils.R +++ b/R/dd_utils.R @@ -516,9 +516,11 @@ sample2 = function(x,size,replace = FALSE,prob = NULL) #' @param trparsopt Initial guess of the parameters to be optimized #' @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})}. @@ -720,6 +722,7 @@ simplex = function(fun, trparsopt, optimpars, ...) #' 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. diff --git a/man/optimizer.Rd b/man/optimizer.Rd index 58553c0..0ef3346 100644 --- a/man/optimizer.Rd +++ b/man/optimizer.Rd @@ -11,7 +11,6 @@ optimizer( fun, trparsopt, jitter = 0, - verbose = FALSE, ... ) } @@ -19,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, @@ -31,13 +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{verbose}{if TRUE, prints intermediate output when using simplex} - -\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 15350c4..ff0f6ef 100644 --- a/man/simplex.Rd +++ b/man/simplex.Rd @@ -4,18 +4,18 @@ \alias{simplex} \title{Carries out optimization using a simplex algorithm (finding a minimum)} \usage{ -simplex(fun, trparsopt, optimpars, verbose = TRUE, ...) +simplex(fun, trparsopt, optimpars, ...) } \arguments{ \item{fun}{Function to be optimized} \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{verbose}{adds verbose output of intermediate evaluations} +\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 optimimized, or settings of the optimization routine} From c592c6310e5e373ef1f07ee22d1e696605f00233 Mon Sep 17 00:00:00 2001 From: Thijs Janzen Date: Tue, 1 Apr 2025 11:19:46 +0200 Subject: [PATCH 12/29] remove duplicate verbose output --- R/dd_utils.R | 14 -------------- 1 file changed, 14 deletions(-) diff --git a/R/dd_utils.R b/R/dd_utils.R index a3e7b76..c8c92c0 100644 --- a/R/dd_utils.R +++ b/R/dd_utils.R @@ -573,13 +573,6 @@ simplex = function(fun, trparsopt, optimpars, ...) cat(string) utils::flush.console() } - - - if (verbose) { - string = paste(string, -fv[1], how, "\n", sep = " ") - cat(string) - utils::flush.console() - } tmp = order(fv) if(numpar == 1) @@ -690,13 +683,6 @@ simplex = function(fun, trparsopt, optimpars, ...) utils::flush.console() } - - if (verbose) { - 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) From a1c6790ff810b4132f035b7f55b764978e97ffe7 Mon Sep 17 00:00:00 2001 From: Thijs Janzen Date: Wed, 2 Apr 2025 11:22:28 +0200 Subject: [PATCH 13/29] Update dd_utils.R --- R/dd_utils.R | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/R/dd_utils.R b/R/dd_utils.R index c8c92c0..986e6ff 100644 --- a/R/dd_utils.R +++ b/R/dd_utils.R @@ -571,7 +571,7 @@ simplex = function(fun, trparsopt, optimpars, ...) } string = paste(string, -fv[1], how, "\n", sep = " ") cat(string) - utils::flush.console() + #utils::flush.console() } tmp = order(fv) @@ -680,7 +680,6 @@ simplex = function(fun, trparsopt, optimpars, ...) } 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)) From 619901b856274c07e9bffd5e8a2842945701fb12 Mon Sep 17 00:00:00 2001 From: rsetienne Date: Fri, 6 Jun 2025 16:18:02 +0200 Subject: [PATCH 14/29] introduce abstol and reltol in pars2 --- R/dd_loglik.R | 4 ++-- R/dd_utils.R | 16 ++++++++-------- 2 files changed, 10 insertions(+), 10 deletions(-) diff --git a/R/dd_loglik.R b/R/dd_loglik.R index 335fa77..473f7fa 100644 --- a/R/dd_loglik.R +++ b/R/dd_loglik.R @@ -181,8 +181,8 @@ 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['abstol'])) abstol <- 1e-10 + if(is.na(pars2['reltol'])) reltol <- 1e-8 brts = -sort(abs(as.numeric(brts)),decreasing = TRUE) if(sum(brts == 0) == 0) { diff --git a/R/dd_utils.R b/R/dd_utils.R index 93429b8..b138ba0 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]) } From 14a70b5f9254fa6d1624990058aa9f0c7b8ec5ee Mon Sep 17 00:00:00 2001 From: rsetienne Date: Fri, 6 Jun 2025 16:44:04 +0200 Subject: [PATCH 15/29] tolint in dd_ML --- R/dd_ML.R | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/R/dd_ML.R b/R/dd_ML.R index e11240b..d81bc52 100644 --- a/R/dd_ML.R +++ b/R/dd_ML.R @@ -145,6 +145,7 @@ dd_ML = function( btorph = 1, soc = 2, tol = c(1E-3, 1E-4, 1E-6), + tolint = c(1E-10,1E-8), maxiter = 1000 * round((1.25)^length(idparsopt)), changeloglikifnoconv = FALSE, optimmethod = 'simplex', @@ -183,7 +184,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,abstol = tolint[1],reltol = tolint[2]) 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") From ee5e4d9fdc48f588f3c5d3316b6483ff26a0f168 Mon Sep 17 00:00:00 2001 From: rsetienne Date: Fri, 6 Jun 2025 16:46:19 +0200 Subject: [PATCH 16/29] Documentation --- R/dd_ML.R | 2 ++ man/dd_ML.Rd | 4 ++++ 2 files changed, 6 insertions(+) diff --git a/R/dd_ML.R b/R/dd_ML.R index d81bc52..0ebcaaa 100644 --- a/R/dd_ML.R +++ b/R/dd_ML.R @@ -95,6 +95,8 @@ parsfixdefault = function(ddmodel,brts,missnumspec,idparsopt) #' = 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 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 diff --git a/man/dd_ML.Rd b/man/dd_ML.Rd index 801eb85..eaa79cc 100644 --- a/man/dd_ML.Rd +++ b/man/dd_ML.Rd @@ -18,6 +18,7 @@ dd_ML( btorph = 1, soc = 2, tol = c(0.001, 1e-04, 1e-06), + tolint = c(1e-10, 1e-08), maxiter = 1000 * round((1.25)^length(idparsopt)), changeloglikifnoconv = FALSE, optimmethod = "simplex", @@ -100,6 +101,9 @@ the phylogeny (1)} 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{maxiter}{Sets the maximum number of iterations in the optimization} \item{changeloglikifnoconv}{if TRUE the loglik will be set to -Inf if ML From 11fbc07cce10280a3df4cf25076fb8f3d0de8caf Mon Sep 17 00:00:00 2001 From: rsetienne Date: Fri, 6 Jun 2025 17:00:28 +0200 Subject: [PATCH 17/29] abstolint and reltolint --- R/dd_ML.R | 2 +- R/dd_loglik.R | 28 ++++++++++++++-------------- 2 files changed, 15 insertions(+), 15 deletions(-) diff --git a/R/dd_ML.R b/R/dd_ML.R index 0ebcaaa..44541d7 100644 --- a/R/dd_ML.R +++ b/R/dd_ML.R @@ -186,7 +186,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,abstol = tolint[1],reltol = tolint[2]) + pars2 = c(res,ddmodel,cond,btorph,verbose,soc,tol,maxiter,abstolint = tolint[1],reltolint = tolint[2]) 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 473f7fa..9bb5101 100644 --- a/R/dd_loglik.R +++ b/R/dd_loglik.R @@ -181,8 +181,8 @@ 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 { - if(is.na(pars2['abstol'])) abstol <- 1e-10 - if(is.na(pars2['reltol'])) reltol <- 1e-8 + if(is.na(pars2['abstolint'])) abstolint <- 1e-10 + if(is.na(pars2['reltolint'])) reltolint <- 1e-8 brts = -sort(abs(as.numeric(brts)),decreasing = TRUE) if(sum(brts == 0) == 0) { @@ -212,7 +212,7 @@ dd_loglik1 = function(pars1,pars2,brts,missnumspec,methode = 'odeint::runge_kutt 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) + 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) { @@ -232,7 +232,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) { @@ -256,7 +256,7 @@ 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 } @@ -271,14 +271,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]) } } @@ -350,8 +350,8 @@ if((ddep == 1) & ((mu == 0 & missnumspec == 0 & floor(K) != ceiling(K) & la > 0. { 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 +if(is.na(pars2['abstolint'])) abstolint <- 1e-16 +if(is.na(pars2['reltolint'])) reltolint <- 1e-10 brts = -sort(abs(as.numeric(brts)),decreasing = TRUE) if(sum(brts == 0) == 0) { @@ -379,7 +379,7 @@ if((mu == 0 & (ddep == 2 | ddep == 2.1 | ddep == 2.2)) | (la == 0 & (ddep == 4 | 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) + #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) @@ -401,7 +401,7 @@ if((mu == 0 & (ddep == 2 | ddep == 2.1 | ddep == 2.2)) | (la == 0 & (ddep == 4 | 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) + #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) @@ -427,7 +427,7 @@ if((mu == 0 & (ddep == 2 | ddep == 2.1 | ddep == 2.2)) | (la == 0 & (ddep == 4 | 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); + #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 } @@ -441,7 +441,7 @@ if((mu == 0 & (ddep == 2 | ddep == 2.1 | ddep == 2.2)) | (la == 0 & (ddep == 4 | #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) + #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 @@ -457,7 +457,7 @@ if((mu == 0 & (ddep == 2 | ddep == 2.1 | ddep == 2.2)) | (la == 0 & (ddep == 4 | #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] From 53f5e6da0129f458c196eb6f8e52ec15124a21fd Mon Sep 17 00:00:00 2001 From: rsetienne Date: Fri, 6 Jun 2025 17:08:00 +0200 Subject: [PATCH 18/29] abstolint, reltolint --- R/dd_loglik.R | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/R/dd_loglik.R b/R/dd_loglik.R index 9bb5101..c94c9fc 100644 --- a/R/dd_loglik.R +++ b/R/dd_loglik.R @@ -181,8 +181,8 @@ 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 { - if(is.na(pars2['abstolint'])) abstolint <- 1e-10 - if(is.na(pars2['reltolint'])) reltolint <- 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'] brts = -sort(abs(as.numeric(brts)),decreasing = TRUE) if(sum(brts == 0) == 0) { @@ -350,8 +350,8 @@ if((ddep == 1) & ((mu == 0 & missnumspec == 0 & floor(K) != ceiling(K) & la > 0. { 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 -if(is.na(pars2['reltolint'])) reltolint <- 1e-10 +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) { From 601f58d1457e6a252233b4a88c61f2d2d68480a1 Mon Sep 17 00:00:00 2001 From: rsetienne Date: Sat, 7 Jun 2025 07:33:23 +0200 Subject: [PATCH 19/29] No checkprobs for normalization --- R/dd_loglik.R | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/R/dd_loglik.R b/R/dd_loglik.R index c94c9fc..08ee548 100644 --- a/R/dd_loglik.R +++ b/R/dd_loglik.R @@ -261,9 +261,11 @@ dd_loglik1 = function(pars1,pars2,brts,missnumspec,methode = 'odeint::runge_kutt 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) { From 128cd6193a5f019ab5cca32fb15fbe0fa868a14d Mon Sep 17 00:00:00 2001 From: rsetienne Date: Sun, 8 Jun 2025 17:16:22 +0200 Subject: [PATCH 20/29] new branch --- DDD.Rproj | 1 - R/RcppExports.R | 7 + R/dd_loglik.R | 379 +++++++++++++++++--------------- R/dd_loglik_rhs.R | 15 ++ src/RcppExports.cpp | 17 ++ src/dd_integrate_log_odeint.cpp | 79 +++++++ 6 files changed, 318 insertions(+), 180 deletions(-) create mode 100644 src/dd_integrate_log_odeint.cpp 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_loglik.R b/R/dd_loglik.R index 08ee548..c2c10f6 100644 --- a/R/dd_loglik.R +++ b/R/dd_loglik.R @@ -132,21 +132,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') @@ -212,8 +212,15 @@ dd_loglik1 = function(pars1,pars2,brts,missnumspec,methode = 'odeint::runge_kutt for(k in 2:(S + 2 - soc)) { k1 = k + (soc - 2) + if(k == 15) { + print(probs[1]) + rhs_func_name <- 'dd_loglik_log_rhs' + probs <- log(probs) + print(probs[1]) + } 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)] + print(c(k,probs[1])) if(is.na(sum(probs)) && pars1[2]/pars1[1] < 1E-4 && missnumspec == 0) { loglik = dd_loglik_high_lambda(pars1 = pars1,pars2 = pars2,brts = brts) @@ -222,9 +229,14 @@ 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(any(flavec(ddep,la,mu,K,r,lx,k1) == 0)) print(which.max(flavec(ddep,la,mu,K,r,lx,k1) == 0)) + if(k >= 15) { + probs <- log(flavec(ddep,la,mu,K,r,lx,k1)) + probs # speciation event + } 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) @@ -238,16 +250,21 @@ dd_loglik1 = function(pars1,pars2,brts,missnumspec,methode = 'odeint::runge_kutt { probs = c(flavec(ddep,la,mu,K,r,lx,k1-1),1) * probs # speciation event } - cp <- check_probs(loglik,probs[1:lx],verbose); loglik <- cp[[1]]; probs[1:lx] <- cp[[2]]; + #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 < 15 & 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 < 15) { + 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 if(cond == 1 | cond == 2) { @@ -256,16 +273,20 @@ dd_loglik1 = function(pars1,pars2,brts,missnumspec,methode = 'odeint::runge_kutt k = soc t1 = brts[1] t2 = brts[S + 2 - soc] + #interval <- seq(t1,t2,(t2 - t1)/1000) + #for(i in 1:999) { + # y <- dd_integrate(probsn,c(interval[i],interval[i + 1]),rhs_func_name,c(pars1,k,ddep),rtol = reltolint,atol = abstolint,method = methode); + # probsn <- y[2,2:(lx+1)] + # cp <- check_probs(logliknorm,probsn,verbose); logliknorm <- cp[[1]]; probsn <- cp[[2]]; + #} 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])} - if(cond == 1) { logliknorm = log(sum(probsc)) } - if(cond == 2) { logliknorm = log(probsc[S + missnumspec - soc + 1])} + if(cond == 1) { logliknorm = logliknorm + log(sum(probsc)) } + if(cond == 2) { logliknorm = logliknorm + log(probsc[S + missnumspec - soc + 1])} } if(cond == 3) { @@ -315,146 +336,147 @@ 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 + print(sum(probsc)) + 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 +490,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 +546,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_bw_rhs' || rhs_func_name == 'dd_loglik_log_rhs') { parsvec = c(dd_loglik_rhs_precomp(pars,initprobs),pars[length(pars) - 1]) } else @@ -551,7 +573,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 ) @@ -597,6 +620,4 @@ dd_ode_odeint = function(initprobs, tvec, parsvec, atol, rtol, method, rhs_name) y[2,2:(lx+1)] <- fun(initprobs, tvec, parsvec, atol, rtol, method) } return(y) -} - - +} \ No newline at end of file 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/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); +} From 426ae18d181b3805e70c63be4501fb0d99608e8b Mon Sep 17 00:00:00 2001 From: rsetienne Date: Mon, 9 Jun 2025 09:25:56 +0200 Subject: [PATCH 21/29] k_threshold --- R/dd_loglik.R | 31 +++++++++++++++++++------------ 1 file changed, 19 insertions(+), 12 deletions(-) diff --git a/R/dd_loglik.R b/R/dd_loglik.R index c2c10f6..c3da97f 100644 --- a/R/dd_loglik.R +++ b/R/dd_loglik.R @@ -151,6 +151,7 @@ dd_loglik = function(pars1,pars2,brts,missnumspec,methode = 'analytical') dd_loglik1 = function(pars1,pars2,brts,missnumspec,methode = 'odeint::runge_kutta_cash_karp54',rhs_func_name = 'dd_loglik_rhs') { + k_threshold <- 10000 if(length(pars2) == 4) { pars2[5] = 0 @@ -212,15 +213,17 @@ dd_loglik1 = function(pars1,pars2,brts,missnumspec,methode = 'odeint::runge_kutt for(k in 2:(S + 2 - soc)) { k1 = k + (soc - 2) - if(k == 15) { - print(probs[1]) - rhs_func_name <- 'dd_loglik_log_rhs' - probs <- log(probs) - print(probs[1]) + 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)] - print(c(k,probs[1])) if(is.na(sum(probs)) && pars1[2]/pars1[1] < 1E-4 && missnumspec == 0) { loglik = dd_loglik_high_lambda(pars1 = pars1,pars2 = pars2,brts = brts) @@ -229,10 +232,14 @@ dd_loglik1 = function(pars1,pars2,brts,missnumspec,methode = 'odeint::runge_kutt } if(k < (S + 2 - soc)) { - if(any(flavec(ddep,la,mu,K,r,lx,k1) == 0)) print(which.max(flavec(ddep,la,mu,K,r,lx,k1) == 0)) - if(k >= 15) { - probs <- log(flavec(ddep,la,mu,K,r,lx,k1)) + probs # speciation event - } else { + 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 } } @@ -253,12 +260,12 @@ 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((k < 15 & 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 { - if(k < 15) { + 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 { From e92b4b8346eb775e4a79a73bbef6058da4ab8346 Mon Sep 17 00:00:00 2001 From: rsetienne Date: Mon, 9 Jun 2025 09:32:02 +0200 Subject: [PATCH 22/29] back to normal? --- R/dd_loglik.R | 36 ++++++++++++++++++------------------ 1 file changed, 18 insertions(+), 18 deletions(-) diff --git a/R/dd_loglik.R b/R/dd_loglik.R index c3da97f..adb8b52 100644 --- a/R/dd_loglik.R +++ b/R/dd_loglik.R @@ -213,15 +213,15 @@ dd_loglik1 = function(pars1,pars2,brts,missnumspec,methode = 'odeint::runge_kutt 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) - } - } + #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) @@ -232,16 +232,16 @@ dd_loglik1 = function(pars1,pars2,brts,missnumspec,methode = 'odeint::runge_kutt } if(k < (S + 2 - soc)) { - 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 - { + #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]]; } From 21e6edc1590aa846eb0b211df0c18938d7e4a0db Mon Sep 17 00:00:00 2001 From: rsetienne Date: Mon, 9 Jun 2025 12:24:40 +0200 Subject: [PATCH 23/29] integration of logprobs --- R/dd_loglik.R | 65 ++++++++++++++++++++++++++------------------------- 1 file changed, 33 insertions(+), 32 deletions(-) diff --git a/R/dd_loglik.R b/R/dd_loglik.R index adb8b52..e810fcb 100644 --- a/R/dd_loglik.R +++ b/R/dd_loglik.R @@ -151,7 +151,7 @@ dd_loglik = function(pars1,pars2,brts,missnumspec,methode = 'analytical') dd_loglik1 = function(pars1,pars2,brts,missnumspec,methode = 'odeint::runge_kutta_cash_karp54',rhs_func_name = 'dd_loglik_rhs') { - k_threshold <- 10000 + k_threshold <- 10 if(length(pars2) == 4) { pars2[5] = 0 @@ -210,18 +210,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) - # } - #} + 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) @@ -232,16 +234,16 @@ dd_loglik1 = function(pars1,pars2,brts,missnumspec,methode = 'odeint::runge_kutt } if(k < (S + 2 - soc)) { - #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 - #{ + 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]]; } @@ -257,7 +259,7 @@ dd_loglik1 = function(pars1,pars2,brts,missnumspec,methode = 'odeint::runge_kutt { probs = c(flavec(ddep,la,mu,K,r,lx,k1-1),1) * probs # speciation event } - #cp <- check_probs(loglik,probs[1:lx],verbose); loglik <- cp[[1]]; probs[1:lx] <- cp[[2]]; + cp <- check_probs(loglik,probs[1:lx],verbose); loglik <- cp[[1]]; probs[1:lx] <- cp[[2]]; } } if((k < k_threshold & probs[1 + missnumspec] <= 0) | loglik == -Inf | is.na(loglik) | is.nan(loglik)) @@ -269,10 +271,12 @@ dd_loglik1 = function(pars1,pars2,brts,missnumspec,methode = 'odeint::runge_kutt 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) + 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) @@ -280,20 +284,16 @@ dd_loglik1 = function(pars1,pars2,brts,missnumspec,methode = 'odeint::runge_kutt k = soc t1 = brts[1] t2 = brts[S + 2 - soc] - #interval <- seq(t1,t2,(t2 - t1)/1000) - #for(i in 1:999) { - # y <- dd_integrate(probsn,c(interval[i],interval[i + 1]),rhs_func_name,c(pars1,k,ddep),rtol = reltolint,atol = abstolint,method = methode); - # probsn <- y[2,2:(lx+1)] - # cp <- check_probs(logliknorm,probsn,verbose); logliknorm <- cp[[1]]; probsn <- cp[[2]]; - #} 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])} + #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) { @@ -463,7 +463,6 @@ dd_loglik2 = function(pars1,pars2,brts,missnumspec) if(soc == 1) { aux = 1:lx } if(soc == 2) { aux = (2:(lx+1)) * (3:(lx+2))/6 } probsc = probsn/aux - print(sum(probsc)) if(cond == 1) { logliknorm = log(sum(probsc)) } if(cond == 2) { logliknorm = log(probsc[S + missnumspec - soc + 1])} } @@ -559,7 +558,7 @@ dd_integrate <- function(initprobs,tvec,rhs_func,pars,rtol,atol,method) rhs_func = match.fun(rhs_func) #} #} - if(rhs_func_name == 'dd_loglik_rhs' || rhs_func_name == 'dd_loglik_bw_rhs' || rhs_func_name == 'dd_loglik_log_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 @@ -627,4 +626,6 @@ dd_ode_odeint = function(initprobs, tvec, parsvec, atol, rtol, method, rhs_name) y[2,2:(lx+1)] <- fun(initprobs, tvec, parsvec, atol, rtol, method) } return(y) -} \ No newline at end of file +} + + From b37758be660a5c0e8c558d248ea39aa8633bb418 Mon Sep 17 00:00:00 2001 From: rsetienne Date: Mon, 9 Jun 2025 13:54:43 +0200 Subject: [PATCH 24/29] Flexibility in log-integration --- R/dd_ML.R | 3 ++- R/dd_loglik.R | 4 ++-- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/R/dd_ML.R b/R/dd_ML.R index 44541d7..9ef9d33 100644 --- a/R/dd_ML.R +++ b/R/dd_ML.R @@ -148,6 +148,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 +187,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 e810fcb..2c8b131 100644 --- a/R/dd_loglik.R +++ b/R/dd_loglik.R @@ -151,7 +151,6 @@ dd_loglik = function(pars1,pars2,brts,missnumspec,methode = 'analytical') dd_loglik1 = function(pars1,pars2,brts,missnumspec,methode = 'odeint::runge_kutta_cash_karp54',rhs_func_name = 'dd_loglik_rhs') { - k_threshold <- 10 if(length(pars2) == 4) { pars2[5] = 0 @@ -184,6 +183,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) { @@ -243,9 +243,9 @@ dd_loglik1 = function(pars1,pars2,brts,missnumspec,methode = 'odeint::runge_kutt } 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) From bc6eada845c988cf0911c4dd57049899d620ca99 Mon Sep 17 00:00:00 2001 From: rsetienne Date: Mon, 9 Jun 2025 21:49:02 +0200 Subject: [PATCH 25/29] Documentation --- R/dd_ML.R | 2 ++ R/dd_loglik.R | 5 +++++ man/dd_ML.Rd | 4 ++++ man/dd_loglik.Rd | 7 ++++++- 4 files changed, 17 insertions(+), 1 deletion(-) diff --git a/R/dd_ML.R b/R/dd_ML.R index 9ef9d33..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 diff --git a/R/dd_loglik.R b/R/dd_loglik.R index 2c8b131..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 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} From 1a2d9863cba05ab1d5d70867ecfda2379e894e6b Mon Sep 17 00:00:00 2001 From: rsetienne Date: Tue, 10 Jun 2025 23:25:29 +0200 Subject: [PATCH 26/29] Avoid error in matrix exponentiation --- R/dd_loglik_M.R | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) 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) } From b49dafc4c033cc06d1b29eb2b79ce0a3cd708597 Mon Sep 17 00:00:00 2001 From: rsetienne Date: Sat, 21 Jun 2025 21:27:41 +0200 Subject: [PATCH 27/29] Change default k_threshold --- R/dd_ML.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/dd_ML.R b/R/dd_ML.R index e9cc9df..4d59d18 100644 --- a/R/dd_ML.R +++ b/R/dd_ML.R @@ -150,7 +150,7 @@ dd_ML = function( soc = 2, tol = c(1E-3, 1E-4, 1E-6), tolint = c(1E-10,1E-8), - k_threshold = 3, + k_threshold = 100, maxiter = 1000 * round((1.25)^length(idparsopt)), changeloglikifnoconv = FALSE, optimmethod = 'simplex', From 29a4d27b224ac44b5fb82019ec2e9c2cb4f7fb66 Mon Sep 17 00:00:00 2001 From: rsetienne Date: Mon, 23 Jun 2025 10:15:09 +0200 Subject: [PATCH 28/29] Documentation update and importing from optim --- R/dd_utils.R | 16 ++++++++-------- man/dd_ML.Rd | 2 +- 2 files changed, 9 insertions(+), 9 deletions(-) diff --git a/R/dd_utils.R b/R/dd_utils.R index 6d31b8a..3577199 100644 --- a/R/dd_utils.R +++ b/R/dd_utils.R @@ -800,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)) diff --git a/man/dd_ML.Rd b/man/dd_ML.Rd index 4da3b84..3ed2c67 100644 --- a/man/dd_ML.Rd +++ b/man/dd_ML.Rd @@ -19,7 +19,7 @@ dd_ML( soc = 2, tol = c(0.001, 1e-04, 1e-06), tolint = c(1e-10, 1e-08), - k_threshold = 3, + k_threshold = 100, maxiter = 1000 * round((1.25)^length(idparsopt)), changeloglikifnoconv = FALSE, optimmethod = "simplex", From dc53d04e8341d12c1df7420399352045b5f13044 Mon Sep 17 00:00:00 2001 From: rsetienne Date: Tue, 24 Jun 2025 12:58:02 +0200 Subject: [PATCH 29/29] Replacing k_threshold by probs_threshold --- R/dd_ML.R | 21 +++++++++++---------- R/dd_loglik.R | 18 +++++++++--------- man/dd_ML.Rd | 19 ++++++++++--------- 3 files changed, 30 insertions(+), 28 deletions(-) diff --git a/R/dd_ML.R b/R/dd_ML.R index 4d59d18..7c1fd4c 100644 --- a/R/dd_ML.R +++ b/R/dd_ML.R @@ -91,14 +91,15 @@ 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 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 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 @@ -150,7 +151,7 @@ dd_ML = function( soc = 2, tol = c(1E-3, 1E-4, 1E-6), tolint = c(1E-10,1E-8), - k_threshold = 100, + probs_threshold = 0, maxiter = 1000 * round((1.25)^length(idparsopt)), changeloglikifnoconv = FALSE, optimmethod = 'simplex', @@ -189,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,abstolint = tolint[1],reltolint = tolint[2],k_threshold = k_threshold) + 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_loglik.R b/R/dd_loglik.R index c1bf988..3fe4318 100644 --- a/R/dd_loglik.R +++ b/R/dd_loglik.R @@ -188,7 +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'] + 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) { @@ -216,17 +216,17 @@ 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) probs <- probs[1:lx] - if(k >= k_threshold) { + if(probs[1] > 0 & probs[1] < probs_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) - } + 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) @@ -239,7 +239,7 @@ dd_loglik1 = function(pars1,pars2,brts,missnumspec,methode = 'odeint::runge_kutt } if(k < (S + 2 - soc)) { - if(k >= k_threshold) { + if(logtrafo) { fac <- flavec(ddep,la,mu,K,r,lx,k1) lx <- min(lx,which(fac == 0) - 1) probs <- probs[1:lx] @@ -267,12 +267,12 @@ 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((k < k_threshold & 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 { - if(k < k_threshold) { + 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 { diff --git a/man/dd_ML.Rd b/man/dd_ML.Rd index 3ed2c67..198514a 100644 --- a/man/dd_ML.Rd +++ b/man/dd_ML.Rd @@ -19,7 +19,7 @@ dd_ML( soc = 2, tol = c(0.001, 1e-04, 1e-06), tolint = c(1e-10, 1e-08), - k_threshold = 100, + probs_threshold = 0, maxiter = 1000 * round((1.25)^length(idparsopt)), changeloglikifnoconv = FALSE, optimmethod = "simplex", @@ -97,16 +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{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{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}