diff --git a/DESCRIPTION b/DESCRIPTION index de0fab8..a006a44 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,8 +1,8 @@ Package: mkin Type: Package Title: Kinetic Evaluation of Chemical Degradation Data -Version: 1.2.0 -Date: 2022-11-16 +Version: 1.3.0 +Date: 2022-11-15 Authors@R: c( person("Johannes", "Ranke", role = c("aut", "cre", "cph"), email = "johannes.ranke@jrwb.de", @@ -22,8 +22,8 @@ Description: Calculation routines based on the FOCUS Kinetics Report (2006, Ranke et al. (2021) . Please note that no warranty is implied for correctness of results or fitness for a particular purpose. -Depends: R (>= 2.15.1), -Imports: stats, graphics, methods, parallel, deSolve, R6, inline (>= 0.3.19), +Depends: R (>= 2.15.1), deSolve +Imports: stats, graphics, methods, parallel, R6, inline (>= 0.3.19), numDeriv, lmtest, pkgbuild, nlme (>= 3.1-151), saemix (>= 3.1), rlang, vctrs Suggests: knitr, rbenchmark, tikzDevice, testthat, rmarkdown, covr, vdiffr, benchmarkme, tibble, stats4, readxl diff --git a/NAMESPACE b/NAMESPACE index 107ffc5..4a41acc 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -147,7 +147,6 @@ export(status) export(tex_listing) export(transform_odeparms) export(which.best) -import(deSolve) import(graphics) import(nlme) importFrom(R6,R6Class) diff --git a/NEWS.md b/NEWS.md index 22d5071..846c7c5 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,11 +1,5 @@ # mkin 1.2.0 (unreleased) -- 'R/saem.R': 'logLik', 'update' and 'anova' methods for 'saem.mmkin' objects. - -- 'R/saem.R': Automatic estimation of start parameters for random effects for the case of mkin transformations, nicely improving convergence and reducing problems with iterative ODE solutions. - -- 'R/status.R': New generic to show status information for fit array objects with methods for 'mmkin', 'mhmkin' and 'multistart' objects. - - 'R/mhmkin.R': New method for performing multiple hierarchical mkin fits in one function call, optionally in parallel. - 'R/mhmkin.R': 'anova.mhmkin' for conveniently comparing the resulting fits. @@ -14,6 +8,12 @@ - 'R/multistart.R': New method for testing multiple start parameters for hierarchical model fits, with function 'llhist' and new generic 'parplot' for diagnostics, and new generics 'which.best' and 'best' for extracting the fit with the highest likelihood +- 'R/saem.R': 'logLik', 'update' and 'anova' methods for 'saem.mmkin' objects. + +- 'R/saem.R': Automatic estimation of start parameters for random effects for the case of mkin transformations, nicely improving convergence and reducing problems with iterative ODE solutions. + +- 'R/status.R': New generic to show status information for fit array objects with methods for 'mmkin', 'mhmkin' and 'multistart' objects. + - 'R/summary.mmkin.R': Summary method for mmkin objects. - 'R/saem.R': Implement and test saemix transformations for FOMC and HS. Also, error out if saemix transformations are requested but not supported. diff --git a/R/mkinfit.R b/R/mkinfit.R index 693778f..0d9246d 100644 --- a/R/mkinfit.R +++ b/R/mkinfit.R @@ -500,6 +500,15 @@ mkinfit <- function(mkinmod, observed, } } + # Get native symbol before iterations info for speed + call_lsoda <- getNativeSymbolInfo("call_lsoda", PACKAGE = "deSolve") + if (solution_type == "deSolve" & use_compiled[1] != FALSE) { + mkinmod$diffs_address <- getNativeSymbolInfo("diffs", + PACKAGE = mkinmod$dll_info[["name"]])$address + mkinmod$initpar_address <- getNativeSymbolInfo("initpar", + PACKAGE = mkinmod$dll_info[["name"]])$address + } + # Get the error model and the algorithm for fitting err_mod <- match.arg(error_model) error_model_algorithm = match.arg(error_model_algorithm) @@ -610,7 +619,8 @@ mkinfit <- function(mkinmod, observed, solution_type = solution_type, use_compiled = use_compiled, method.ode = method.ode, - atol = atol, rtol = rtol, ...) + atol = atol, rtol = rtol, + call_lsoda = call_lsoda, ...) observed_index <- cbind(as.character(observed$time), as.character(observed$name)) observed$predicted <- out[observed_index] @@ -892,7 +902,10 @@ mkinfit <- function(mkinmod, observed, fit$calls <- calls fit$time <- fit_time - # We also need the model and a model name for summary and plotting + # We also need the model and a model name for summary and plotting, + # but without address info that will become invalid + mkinmod$diffs_address <- NULL + mkinmod$initpar_address <- NULL fit$mkinmod <- mkinmod fit$mkinmod$name <- mkinmod_name fit$obs_vars <- obs_vars diff --git a/R/mkinmod.R b/R/mkinmod.R index d8740ae..b1fb57c 100644 --- a/R/mkinmod.R +++ b/R/mkinmod.R @@ -45,7 +45,9 @@ #' @param dll_dir Directory where an DLL object, if generated internally by #' [inline::cfunction()], should be saved. The DLL will only be stored in a #' permanent location for use in future sessions, if 'dll_dir' and 'name' -#' are specified. +#' are specified. This is helpful if fit objects are cached e.g. by knitr, +#' as the cache remains functional across sessions if the DLL is stored in +#' a user defined location. #' @param unload If a DLL from the target location in 'dll_dir' is already #' loaded, should that be unloaded first? #' @param overwrite If a file exists at the target DLL location in 'dll_dir', diff --git a/R/mkinpredict.R b/R/mkinpredict.R index 0dc9cf5..11a3b35 100644 --- a/R/mkinpredict.R +++ b/R/mkinpredict.R @@ -19,26 +19,24 @@ #' @param solution_type The method that should be used for producing the #' predictions. This should generally be "analytical" if there is only one #' observed variable, and usually "deSolve" in the case of several observed -#' variables. The third possibility "eigen" is faster but not applicable to -#' some models e.g. using FOMC for the parent compound. +#' variables. The third possibility "eigen" is fast in comparison to uncompiled +#' ODE models, but not applicable to some models, e.g. using FOMC for the +#' parent compound. #' @param method.ode The solution method passed via [mkinpredict] to [ode]] in -#' case the solution type is "deSolve". The default "lsoda" is performant, but -#' sometimes fails to converge. +#' case the solution type is "deSolve" and we are not using compiled code. #' @param use_compiled If set to \code{FALSE}, no compiled version of the #' [mkinmod] model is used, even if is present. -#' @param atol Absolute error tolerance, passed to [ode]. Default is 1e-8, -#' lower than in [lsoda]. -#' @param rtol Absolute error tolerance, passed to [ode]. Default is 1e-10, -#' much lower than in [lsoda]. -#' @param maxsteps Maximum number of steps, passed to [ode]. +#' @param atol Absolute error tolerance, passed to the ode solver. +#' @param rtol Absolute error tolerance, passed to the ode solver. +#' @param maxsteps Maximum number of steps, passed to the ode solver. #' @param map_output Boolean to specify if the output should list values for #' the observed variables (default) or for all state variables (if set to #' FALSE). Setting this to FALSE has no effect for analytical solutions, #' as these always return mapped output. #' @param na_stop Should it be an error if [ode] returns NaN values +#' @param call_lsoda The address of the compiled function "call_lsoda" #' @param \dots Further arguments passed to the ode solver in case such a #' solver is used. -#' @import deSolve #' @return A matrix with the numeric solution in wide format #' @author Johannes Ranke #' @examples @@ -116,9 +114,10 @@ mkinpredict.mkinmod <- function(x, outtimes = seq(0, 120, by = 0.1), solution_type = "deSolve", use_compiled = "auto", - method.ode = "lsoda", atol = 1e-8, rtol = 1e-10, maxsteps = 20000, + method.ode = "lsoda", atol = 1e-8, rtol = 1e-10, maxsteps = 20000L, map_output = TRUE, na_stop = TRUE, + call_lsoda = NULL, ...) { @@ -173,20 +172,80 @@ mkinpredict.mkinmod <- function(x, if (solution_type == "deSolve") { if (!is.null(x$cf) & use_compiled[1] != FALSE) { - out <- deSolve::ode( - y = odeini, - times = outtimes, - func = "diffs", - initfunc = "initpar", - dllname = if (is.null(x$dll_info)) inline::getDynLib(x$cf)[["name"]] - else x$dll_info[["name"]], - parms = odeparms[x$parms], # Order matters when using compiled models - method = method.ode, - atol = atol, - rtol = rtol, - maxsteps = maxsteps, - ... +# out <- deSolve::ode( +# y = odeini, +# times = outtimes, +# func = "diffs", +# initfunc = "initpar", +# dllname = x$dll_info[["name"]], +# parms = odeparms[x$parms], # Order matters when using compiled models +# method = method.ode, +# atol = atol, +# rtol = rtol, +# maxsteps = maxsteps, +# ... +# ) +# + # Prepare call to "call_lsoda" + # Simplified code from deSolve::lsoda() adapted to our use case + if (is.null(call_lsoda)) { + call_lsoda <- getNativeSymbolInfo("call_lsoda", PACKAGE = "deSolve") + } + if (is.null(x$diffs_address)) { + x$diffs_address <- getNativeSymbolInfo("diffs", + PACKAGE = x$dll_info[["name"]])$address + x$initpar_address <- getNativeSymbolInfo("initpar", + PACKAGE = x$dll_info[["name"]])$address + } + rwork <- vector("double", 20) + rwork[] <- 0. + rwork[6] <- max(abs(diff(outtimes))) + + iwork <- vector("integer", 20) + iwork[] <- 0 + iwork[6] <- maxsteps + + n <- length(odeini) + lmat <- n^2 + 2 # from deSolve::lsoda(), for Jacobian type full, internal (2) + # hard-coded default values of lsoda(): + maxordn <- 12L + maxords <- 5L + lrn <- 20 + n * (maxordn + 1) + 3 * n # length in case non-stiff method + lrs <- 20 + n * (maxords + 1) + 3 * n + lmat # length in case stiff method + lrw <- max(lrn, lrs) # actual length: max of both + liw <- 20 + n + + on.exit(.C("unlock_solver", PACKAGE = "deSolve")) + + out_raw <- .Call(call_lsoda, + as.double(odeini), # y + as.double(outtimes), # times + x$diffs_address, # derivfunc + as.double(odeparms[x$parms]), # parms + rtol, atol, + NULL, NULL, # rho, tcrit + NULL, # jacfunc + x$initpar_address, # initfunc + NULL, # eventfunc + 0L, # verbose + 1L, # iTask + as.double(rwork), # rWork + as.integer(iwork), # iWork + 2L, # jT full Jacobian calculated internally + 0L, # nOut + as.integer(lrw), # lRw + as.integer(liw), # lIw + 1L, # Solver + NULL, # rootfunc + 0L, as.double(0), 0L, # nRoot, Rpar, Ipar + 0L, # Type + list(fmat = 0L, tmat = 0L, imat = 0L, ModelForc = NULL), # flist + list(), # elist + list(islag = 0L) # elag ) + + out <- t(out_raw) + colnames(out) <- c("time", mod_vars) } else { mkindiff <- function(t, state, parms) { diff --git a/R/nlme.mmkin.R b/R/nlme.mmkin.R index 09cb84b..e193e5e 100644 --- a/R/nlme.mmkin.R +++ b/R/nlme.mmkin.R @@ -186,10 +186,24 @@ nlme.mmkin <- function(model, data = "auto", thisCall[["control"]] <- control } + # Provide the address of call_lsoda to the fitting function + call_lsoda <- getNativeSymbolInfo("call_lsoda", PACKAGE = "deSolve") + if (model[[1]]$solution_type == "deSolve" & !is.null(model[[1]]$mkinmod$cf)) { + # The mkinmod stored in the first fit will be used by nlme + model[[1]]$mkinmod$diffs_address <- getNativeSymbolInfo("diffs", + PACKAGE = model[[1]]$mkinmod$dll_info[["name"]])$address + model[[1]]$mkinmod$initpar_address <- getNativeSymbolInfo("initpar", + PACKAGE = model[[1]]$mkinmod$dll_info[["name"]])$address + } + fit_time <- system.time(val <- do.call("nlme.formula", thisCall)) val$time <- fit_time val$mkinmod <- model[[1]]$mkinmod + # Don't return addresses that will become invalid + val$mkinmod$diffs_address <- NULL + val$mkinmod$initpar_address <- NULL + val$data <- thisCall[["data"]] val$mmkin <- model if (is.list(start)) val$mean_dp_start <- start$fixed diff --git a/R/read_spreadsheet.R b/R/read_spreadsheet.R index 7ad09c3..a20af6d 100644 --- a/R/read_spreadsheet.R +++ b/R/read_spreadsheet.R @@ -37,7 +37,7 @@ #' and moisture normalisation factors in the sheet 'Datasets'? #' @export read_spreadsheet <- function(path, valid_datasets = "all", - parent_only = FALSE, normalize = TRUE) + parent_only = TRUE, normalize = TRUE) { if (!requireNamespace("readxl", quietly = TRUE)) stop("Please install the readxl package to use this function") diff --git a/R/saem.R b/R/saem.R index 696ea0e..5b8021d 100644 --- a/R/saem.R +++ b/R/saem.R @@ -120,12 +120,12 @@ utils::globalVariables(c("predicted", "std")) #' summary(f_saem_dfop_sfo, data = TRUE) #' #' # The following takes about 6 minutes -#' #f_saem_dfop_sfo_deSolve <- saem(f_mmkin["DFOP-SFO", ], solution_type = "deSolve", -#' # control = list(nbiter.saemix = c(200, 80), nbdisplay = 10)) +#' f_saem_dfop_sfo_deSolve <- saem(f_mmkin["DFOP-SFO", ], solution_type = "deSolve", +#' nbiter.saemix = c(200, 80)) #' -#' #saemix::compare.saemix(list( -#' # f_saem_dfop_sfo$so, -#' # f_saem_dfop_sfo_deSolve$so)) +#' #anova( +#' # f_saem_dfop_sfo, +#' # f_saem_dfop_sfo_deSolve)) #' #' # If the model supports it, we can also use eigenvalue based solutions, which #' # take a similar amount of time @@ -580,6 +580,15 @@ saemix_model <- function(object, solution_type = "auto", transform_rates <- object[[1]]$transform_rates transform_fractions <- object[[1]]$transform_fractions + # Get native symbol info for speed + call_lsoda <- getNativeSymbolInfo("call_lsoda", PACKAGE = "deSolve") + if (solution_type == "deSolve" & !is.null(mkin_model$cf)) { + mkin_model$diffs_address <- getNativeSymbolInfo("diffs", + PACKAGE = mkin_model$dll_info[["name"]])$address + mkin_model$initpar_address <- getNativeSymbolInfo("initpar", + PACKAGE = mkin_model$dll_info[["name"]])$address + } + # Define the model function model_function <- function(psi, id, xidep) { @@ -613,7 +622,8 @@ saemix_model <- function(object, solution_type = "auto", odeparms = odeparms, odeini = odeini, solution_type = solution_type, outtimes = sort(unique(i_time)), - na_stop = FALSE + na_stop = FALSE, + call_lsoda = call_lsoda ) out_index <- cbind(as.character(i_time), as.character(i_name)) diff --git a/docs/dev/articles/mkin.html b/docs/dev/articles/mkin.html index 27e532a..6bfb63b 100644 --- a/docs/dev/articles/mkin.html +++ b/docs/dev/articles/mkin.html @@ -34,7 +34,7 @@ mkin - 1.2.0 + 1.1.0 @@ -44,7 +44,7 @@ Functions and data
  • - Example evaluation of FOCUS Laboratory Data L1 to L3
  • - Example evaluations of dimethenamid data from 2018 with nonlinear mixed-effects models -
    • -
  • - Short demo of the multistart method + Example evaluation of FOCUS Example Dataset Z
  • Performance benefit by using compiled model definitions in mkin
    • -
  • - Example evaluation of FOCUS Example Dataset Z -
  • Calculation of time weighted average concentrations with mkin
    • @@ -78,10 +72,7 @@ Example evaluation of NAFTA SOP Attachment examples
  • - Benchmark timings for mkin -
    • -
  • - Benchmark timings for saem.mmkin + Some benchmark timings
    • @@ -112,7 +103,7 @@

    Introduction to mkin

    Johannes Ranke

    -

    Last change 15 February 2021 (rebuilt 2022-11-16)

    +

    Last change 15 February 2021 (rebuilt 2022-02-28)

    Source: vignettes/mkin.rmd @@ -127,34 +118,34 @@

    In the regulatory evaluation of chemical substances like plant protection products (pesticides), biocides and other chemicals, degradation data play an important role. For the evaluation of pesticide degradation experiments, detailed guidance has been developed, based on nonlinear optimisation. The R add-on package mkin implements fitting some of the models recommended in this guidance from within R and calculates some statistical measures for data series within one or more compartments, for parent and metabolites.

    -library("mkin", quietly = TRUE)
-# Define the kinetic model
-m_SFO_SFO_SFO <- mkinmod(parent = mkinsub("SFO", "M1"),
-                         M1 = mkinsub("SFO", "M2"),
-                         M2 = mkinsub("SFO"),
-                         use_of_ff = "max", quiet = TRUE)
-
-
-# Produce model predictions using some arbitrary parameters
-sampling_times = c(0, 1, 3, 7, 14, 28, 60, 90, 120)
-d_SFO_SFO_SFO <- mkinpredict(m_SFO_SFO_SFO,
-  c(k_parent = 0.03,
-    f_parent_to_M1 = 0.5, k_M1 = log(2)/100,
-    f_M1_to_M2 = 0.9, k_M2 = log(2)/50),
-  c(parent = 100, M1 = 0, M2 = 0),
-  sampling_times)
-
-# Generate a dataset by adding normally distributed errors with
-# standard deviation 3, for two replicates at each sampling time
-d_SFO_SFO_SFO_err <- add_err(d_SFO_SFO_SFO, reps = 2,
-                             sdfunc = function(x) 3,
-                             n = 1, seed = 123456789 )
-
-# Fit the model to the dataset
-f_SFO_SFO_SFO <- mkinfit(m_SFO_SFO_SFO, d_SFO_SFO_SFO_err[[1]], quiet = TRUE)
-
-# Plot the results separately for parent and metabolites
-plot_sep(f_SFO_SFO_SFO, lpos = c("topright", "bottomright", "bottomright"))
    +library("mkin", quietly = TRUE) +# Define the kinetic model +m_SFO_SFO_SFO <- mkinmod(parent = mkinsub("SFO", "M1"), + M1 = mkinsub("SFO", "M2"), + M2 = mkinsub("SFO"), + use_of_ff = "max", quiet = TRUE) + + +# Produce model predictions using some arbitrary parameters +sampling_times = c(0, 1, 3, 7, 14, 28, 60, 90, 120) +d_SFO_SFO_SFO <- mkinpredict(m_SFO_SFO_SFO, + c(k_parent = 0.03, + f_parent_to_M1 = 0.5, k_M1 = log(2)/100, + f_M1_to_M2 = 0.9, k_M2 = log(2)/50), + c(parent = 100, M1 = 0, M2 = 0), + sampling_times) + +# Generate a dataset by adding normally distributed errors with +# standard deviation 3, for two replicates at each sampling time +d_SFO_SFO_SFO_err <- add_err(d_SFO_SFO_SFO, reps = 2, + sdfunc = function(x) 3, + n = 1, seed = 123456789 ) + +# Fit the model to the dataset +f_SFO_SFO_SFO <- mkinfit(m_SFO_SFO_SFO, d_SFO_SFO_SFO_err[[1]], quiet = TRUE) + +# Plot the results separately for parent and metabolites +plot_sep(f_SFO_SFO_SFO, lpos = c("topright", "bottomright", "bottomright"))

    @@ -233,10 +224,10 @@

    Ranke, J. 2021. ‘mkin‘: Kinetic Evaluation of Chemical Degradation Data. https://CRAN.R-project.org/package=mkin.

    -

    Ranke, J., and R. Lehmann. 2012. “Parameter Reliability in Kinetic Evaluation of Environmental Metabolism Data - Assessment and the Influence of Model Specification.” In SETAC World 20-24 May. Berlin. https://jrwb.de/posters/Poster_SETAC_2012_Kinetic_parameter_uncertainty_model_parameterization_Lehmann_Ranke.pdf.

    +

    Ranke, J., and R. Lehmann. 2012. “Parameter Reliability in Kinetic Evaluation of Environmental Metabolism Data - Assessment and the Influence of Model Specification.” In SETAC World 20-24 May. Berlin.

    -

    ———. 2015. “To T-Test or Not to T-Test, That Is the Question.” In XV Symposium on Pesticide Chemistry 2-4 September 2015. Piacenza. https://jrwb.de/posters/piacenza_2015.pdf.

    +

    ———. 2015. “To T-Test or Not to T-Test, That Is the Question.” In XV Symposium on Pesticide Chemistry 2-4 September 2015. Piacenza. http://chem.uft.uni-bremen.de/ranke/posters/piacenza_2015.pdf.

    Ranke, Johannes, and Stefan Meinecke. 2019. “Error Models for the Kinetic Evaluation of Chemical Degradation Data.” Environments 6 (12). https://doi.org/10.3390/environments6120124.

    @@ -271,7 +262,7 @@

    -

    Site built with pkgdown 2.0.6.

    +

    Site built with pkgdown 2.0.2.

    diff --git a/docs/dev/articles/web_only/saem_benchmarks.html b/docs/dev/articles/web_only/saem_benchmarks.html index afff038..e54bc38 100644 --- a/docs/dev/articles/web_only/saem_benchmarks.html +++ b/docs/dev/articles/web_only/saem_benchmarks.html @@ -112,7 +112,7 @@

    Benchmark timings for saem.mmkin

    Johannes Ranke

    -

    Last change 14 November 2022 (rebuilt 2022-11-16)

    +

    Last change 14 November 2022 (rebuilt 2022-11-15)

    Source: vignettes/web_only/saem_benchmarks.rmd @@ -309,10 +309,10 @@ Linux 1.2.0 3.2 -2.156 -4.647 -4.296 -4.951 +2.11 +4.632 +4.264 +4.93

    Two-component error fits for SFO, DFOP, SFORB and HS.

    @@ -332,10 +332,10 @@ Linux 1.2.0 3.2 -5.645 -7.415 -7.848 -7.967 +5.602 +7.373 +7.815 +7.831
    @@ -357,8 +357,8 @@ Linux 1.2.0 3.2 -24.182 -783.932 +24.014 +749.699 @@ -379,7 +379,7 @@ Linux 1.2.0 3.2 -1322.5 +1249.834 diff --git a/docs/dev/news/index.html b/docs/dev/news/index.html index 3353922..8448ebc 100644 --- a/docs/dev/news/index.html +++ b/docs/dev/news/index.html @@ -89,13 +89,13 @@

    mkin 1.2.0 (unreleased)

    -