From efab37957381919c21d874906ce870f4941c760a Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Sat, 9 May 2020 21:18:42 +0200 Subject: Avoid the call to merge for analytical solutions This increases performance up to a factor of five! --- R/create_deg_func.R | 32 ++++++++--------- R/mkinerrmin.R | 28 +++++++-------- R/mkinfit.R | 75 ++++++++++++++++++++-------------------- R/mkinmod.R | 2 +- R/mkinpredict.R | 17 +++++++-- R/plot.mkinfit.R | 1 + man/mkinfit.Rd | 6 ---- man/mkinmod.Rd | 2 +- man/mkinpredict.Rd | 35 +++++++++++-------- man/plot.mkinfit.Rd | 1 + test.log | 33 +++++++++--------- tests/testthat/FOCUS_2006_D.csf | 2 +- tests/testthat/test_analytical.R | 11 ++++++ tests/testthat/test_tests.R | 2 +- 14 files changed, 135 insertions(+), 112 deletions(-) create mode 100644 tests/testthat/test_analytical.R diff --git a/R/create_deg_func.R b/R/create_deg_func.R index fb419a98..886c5e8b 100644 --- a/R/create_deg_func.R +++ b/R/create_deg_func.R @@ -21,19 +21,17 @@ create_deg_func <- function(spec, use_of_ff = c("min", "max")) { use_of_ff <- match.arg(use_of_ff) - min_ff <- use_of_ff == "min" - obs_vars <- names(spec) - n <- character(0) - parent <- obs_vars[1] parent_type <- spec[[1]]$type supported <- TRUE # This may be modified below - n[1] <- paste0(parent, " = ", parent_type, ".solution(outtimes, odeini['", parent, + predicted_text <- character(0) + + predicted_text[parent] <- paste0(parent_type, ".solution(t, odeini['", parent, if (parent_type == "SFORB") "_free", "'], ", switch(parent_type, SFO = paste0("k_", parent, if (min_ff) "_sink" else "", ")"), @@ -59,25 +57,27 @@ create_deg_func <- function(spec, use_of_ff = c("min", "max")) { k2 <- paste0("k_", n2) f12 <- paste0("f_", n1, "_to_", n2) if (parent_type == "SFO") { - n[2] <- paste0(n2, " = (((", k2, "-", k1, ")*", n20, "-", f12, "*", k1, "*", n10, ")*exp(-", k2, "*outtimes)+", - f12, "*", k1, "*", n10, "*exp(-", k1, "*outtimes))/(", k2, "-", k1, ")") + predicted_text[n2] <- paste0( + "(((", k2, "-", k1, ")*", n20, "-", f12, "*", k1, "*", n10, ")*exp(-", k2, "*t)+", + f12, "*", k1, "*", n10, "*exp(-", k1, "*t))/(", k2, "-", k1, ")") } } } if (supported) { - all_n <- paste(n, collapse = ",\n") + deg_func <- function(observed, odeini, odeparms) {} f_body <- paste0("{\n", - "out <- with(as.list(odeparms), {\n", - "data.frame(\n", - "time = outtimes,\n", - all_n, ")\n", + "predicted <- numeric(0)\n", + "with(as.list(odeparms), {\n") + for (obs_var in obs_vars) { + f_body <- paste0(f_body, + "t <- observed[observed$name == '", obs_var, "', 'time']\n", + "predicted <<- c(predicted, ", predicted_text[obs_var], ")\n") + } + f_body <- paste0(f_body, "})\n", - "return(out)\n}\n" - ) - - deg_func <- function(odeini, odeparms, outtimes) {} + "return(predicted)\n}\n") body(deg_func) <- parse(text = f_body) return(deg_func) diff --git a/R/mkinerrmin.R b/R/mkinerrmin.R index 0b647b81..1994aed0 100644 --- a/R/mkinerrmin.R +++ b/R/mkinerrmin.R @@ -1,12 +1,12 @@ if(getRversion() >= '2.15.1') utils::globalVariables(c("name", "value_mean")) #' Calculate the minimum error to assume in order to pass the variance test -#' +#' #' This function finds the smallest relative error still resulting in passing #' the chi-squared test as defined in the FOCUS kinetics report from 2006. -#' +#' #' This function is used internally by \code{\link{summary.mkinfit}}. -#' +#' #' @param fit an object of class \code{\link{mkinfit}}. #' @param alpha The confidence level chosen for the chi-squared test. #' @importFrom stats qchisq aggregate @@ -25,43 +25,41 @@ if(getRversion() >= '2.15.1') utils::globalVariables(c("name", "value_mean")) #' \url{http://esdac.jrc.ec.europa.eu/projects/degradation-kinetics} #' @keywords manip #' @examples -#' +#' #' SFO_SFO = mkinmod(parent = mkinsub("SFO", to = "m1"), #' m1 = mkinsub("SFO"), #' use_of_ff = "max") -#' +#' #' fit_FOCUS_D = mkinfit(SFO_SFO, FOCUS_2006_D, quiet = TRUE) #' round(mkinerrmin(fit_FOCUS_D), 4) #' \dontrun{ #' fit_FOCUS_E = mkinfit(SFO_SFO, FOCUS_2006_E, quiet = TRUE) #' round(mkinerrmin(fit_FOCUS_E), 4) #' } -#' +#' #' @export mkinerrmin <- function(fit, alpha = 0.05) { parms.optim <- fit$par kinerrmin <- function(errdata, n.parms) { - means.mean <- mean(errdata$value_mean, na.rm = TRUE) - df = length(errdata$value_mean) - n.parms + means.mean <- mean(errdata$observed, na.rm = TRUE) + df = nrow(errdata) - n.parms err.min <- sqrt((1 / qchisq(1 - alpha, df)) * - sum((errdata$value_mean - errdata$value_pred)^2)/(means.mean^2)) + sum((errdata$observed - errdata$predicted)^2)/(means.mean^2)) return(list(err.min = err.min, n.optim = n.parms, df = df)) } - means <- aggregate(value ~ time + name, data = fit$observed, mean, na.rm=TRUE) - errdata <- merge(means, fit$predicted, by = c("time", "name"), - suffixes = c("_mean", "_pred")) - errdata <- errdata[order(errdata$time, errdata$name), ] + errdata <- aggregate(cbind(observed, predicted) ~ time + variable, data = fit$data, mean, na.rm=TRUE) + errdata <- errdata[order(errdata$time, errdata$variable), ] # Remove values at time zero for variables whose value for state.ini is fixed, # as these will not have any effect in the optimization and should therefore not # be counted as degrees of freedom. fixed_initials = gsub("_0$", "", rownames(subset(fit$fixed, type = "state"))) - errdata <- subset(errdata, !(time == 0 & name %in% fixed_initials)) + errdata <- subset(errdata, !(time == 0 & variable %in% fixed_initials)) n.optim.overall <- length(parms.optim) - length(fit$errparms) @@ -73,7 +71,7 @@ mkinerrmin <- function(fit, alpha = 0.05) # The degrees of freedom are counted according to FOCUS kinetics (2011, p. 164) for (obs_var in fit$obs_vars) { - errdata.var <- subset(errdata, name == obs_var) + errdata.var <- subset(errdata, variable == obs_var) # Check if initial value is optimised n.initials.optim <- length(grep(paste(obs_var, ".*", "_0", sep=""), names(parms.optim))) diff --git a/R/mkinfit.R b/R/mkinfit.R index 6cfb88c0..c33afbf0 100644 --- a/R/mkinfit.R +++ b/R/mkinfit.R @@ -288,6 +288,10 @@ mkinfit <- function(mkinmod, observed, observed <- subset(observed, value != 0) } + # Sort observed values for efficient analytical predictions + observed$name <- ordered(observed$name, levels = obs_vars) + observed <- observed[order(observed$name, observed$time), ] + # Obtain data for decline from maximum mean value if requested if (from_max_mean) { # This is only used for simple decline models @@ -520,7 +524,7 @@ mkinfit <- function(mkinmod, observed, } # Unique outtimes for model solution. - outtimes = sort(unique(observed$time)) + outtimes <- sort(unique(observed$time)) # Define the objective function for optimisation, including (back)transformations cost_function <- function(P, trans = TRUE, OLS = FALSE, fixed_degparms = FALSE, fixed_errparms = FALSE, update_data = TRUE, ...) @@ -578,47 +582,51 @@ mkinfit <- function(mkinmod, observed, } # Solve the system with current parameter values - out <- mkinpredict(mkinmod, parms, - odeini, outtimes, - solution_type = solution_type, - use_compiled = use_compiled, - method.ode = method.ode, - atol = atol, rtol = rtol, ...) - - out_long <- mkin_wide_to_long(out, time = "time") - - # Surprisingly, the next line is the one taking the most time for one, two - # or three observed variables if we use compiled ODE models - # as evidenced by use of github:hadley/lineprof - cost_data <- merge(observed[c("name", "time", "value")], out_long, - by = c("name", "time"), suffixes = c(".observed", ".predicted")) + if (solution_type == "analytical") { + observed$predicted <- mkinmod$deg_func(observed, odeini, parms) + } else { + out <- mkinpredict(mkinmod, parms, + odeini, outtimes, + solution_type = solution_type, + use_compiled = use_compiled, + method.ode = method.ode, + atol = atol, rtol = rtol, ...) + + out_long <- mkin_wide_to_long(out, time = "time") + + out_merged <- merge(observed[c("name", "time")], out_long) + out_merged$name <- ordered(out_merged$name, levels = obs_vars) + out_merged <- out_merged[order(out_merged$name, out_merged$time), ] + observed$predicted <- out_merged$value + } + # Define standard deviation for each observation if (err_mod == "const") { - cost_data$std <- if (OLS) NA else cost_errparms["sigma"] + observed$std <- if (OLS) NA else cost_errparms["sigma"] } if (err_mod == "obs") { - std_names <- paste0("sigma_", cost_data$name) - cost_data$std <- cost_errparms[std_names] + std_names <- paste0("sigma_", observed$name) + observed$std <- cost_errparms[std_names] } if (err_mod == "tc") { - cost_data$std <- sqrt(cost_errparms["sigma_low"]^2 + cost_data$value.predicted^2 * cost_errparms["rsd_high"]^2) + observed$std <- sqrt(cost_errparms["sigma_low"]^2 + observed$predicted^2 * cost_errparms["rsd_high"]^2) } + # Calculate model cost if (OLS) { # Cost is the sum of squared residuals - cost <- with(cost_data, sum((value.observed - value.predicted)^2)) + cost <- with(observed, sum((value - predicted)^2)) } else { # Cost is the negative log-likelihood - cost <- - with(cost_data, - sum(dnorm(x = value.observed, mean = value.predicted, sd = std, log = TRUE))) + cost <- - with(observed, + sum(dnorm(x = value, mean = predicted, sd = std, log = TRUE))) } # We update the current cost and data during the optimisation, not # during hessian calculations if (update_data) { - assign("out_predicted", out_long, inherits = TRUE) - assign("current_data", cost_data, inherits = TRUE) + assign("current_data", observed, inherits = TRUE) if (cost < cost.current) { assign("cost.current", cost, inherits = TRUE) @@ -718,7 +726,7 @@ mkinfit <- function(mkinmod, observed, degparms <- fit$par # Get the maximum likelihood estimate for sigma at the optimum parameter values - current_data$residual <- current_data$value.observed - current_data$value.predicted + current_data$residual <- current_data$value - current_data$predicted sigma_mle <- sqrt(sum(current_data$residual^2)/nrow(current_data)) # Use that estimate for the constant variance, or as first guess if err_mod = "obs" @@ -855,11 +863,7 @@ mkinfit <- function(mkinmod, observed, # We also need the model and a model name for summary and plotting fit$mkinmod <- mkinmod fit$mkinmod$name <- mkinmod_name - - # We need data and predictions for summary and plotting - fit$observed <- observed fit$obs_vars <- obs_vars - fit$predicted <- out_predicted # Residual sum of squares as a function of the fitted parameters fit$rss <- function(P) cost_function(P, OLS = TRUE, update_data = FALSE) @@ -886,15 +890,10 @@ mkinfit <- function(mkinmod, observed, fit$fixed$type = c(rep("state", length(state.ini.fixed)), rep("deparm", length(parms.fixed))) - # Sort observed, predicted and residuals - current_data$name <- ordered(current_data$name, levels = obs_vars) - - ordered_data <- current_data[order(current_data$name, current_data$time), ] - - fit$data <- data.frame(time = ordered_data$time, - variable = ordered_data$name, - observed = ordered_data$value.observed, - predicted = ordered_data$value.predicted) + fit$data <- data.frame(time = current_data$time, + variable = current_data$name, + observed = current_data$value, + predicted = current_data$predicted) fit$data$residual <- fit$data$observed - fit$data$predicted diff --git a/R/mkinmod.R b/R/mkinmod.R index 21551861..7e9e9248 100644 --- a/R/mkinmod.R +++ b/R/mkinmod.R @@ -70,7 +70,7 @@ #' @examples #' #' # Specify the SFO model (this is not needed any more, as we can now mkinfit("SFO", ...) -#' SFO <- mkinmod(parent = list(type = "SFO")) +#' SFO <- mkinmod(parent = mkinsub("SFO")) #' #' # One parent compound, one metabolite, both single first order #' SFO_SFO <- mkinmod( diff --git a/R/mkinpredict.R b/R/mkinpredict.R index 8c76ca05..75582fac 100644 --- a/R/mkinpredict.R +++ b/R/mkinpredict.R @@ -92,10 +92,15 @@ #' c(parent = 100, m1 = 0), seq(0, 20, by = 0.1), #' solution_type = "analytical", use_compiled = FALSE)[201,]) #' } +#' analytical = mkinpredict(SFO_SFO, +#' c(k_parent = 0.15, f_parent_to_m1 = 0.5, k_m1 = 0.01), +#' c(parent = 100, m1 = 0), seq(0, 20, by = 0.1), +#' solution_type = "analytical", use_compiled = FALSE)[201,] #' #' \dontrun{ #' # Predict from a fitted model #' f <- mkinfit(SFO_SFO, FOCUS_2006_C, quiet = TRUE) +#' f <- mkinfit(SFO_SFO, FOCUS_2006_C, quiet = TRUE, solution_type = "analytical") #' head(mkinpredict(f)) #' } #' @@ -131,8 +136,16 @@ mkinpredict.mkinmod <- function(x, } if (solution_type == "analytical") { - out <- x$deg_func(odeini = odeini, - odeparms = odeparms, outtimes = outtimes) + # This is clumsy, as we wanted fast analytical predictions for mkinfit + out <- data.frame(time = outtimes) + obs_vars <- names(x$spec) + pseudo_observed <- + data.frame(name = rep(obs_vars, each = length(outtimes)), + time = rep(outtimes, length(obs_vars))) + pseudo_observed$predicted <- x$deg_func(pseudo_observed, odeini, odeparms) + for (obs_var in obs_vars) { + out[obs_var] <- pseudo_observed[pseudo_observed$name == obs_var, "predicted"] + } } if (solution_type == "eigen") { diff --git a/R/plot.mkinfit.R b/R/plot.mkinfit.R index 4b5bfa3a..b2cb4890 100644 --- a/R/plot.mkinfit.R +++ b/R/plot.mkinfit.R @@ -69,6 +69,7 @@ if(getRversion() >= '2.15.1') utils::globalVariables(c("type", "variable", "obse #' \dontrun{ #' SFO_SFO <- mkinmod(parent = mkinsub("SFO", "m1", full = "Parent"), #' m1 = mkinsub("SFO", full = "Metabolite M1" )) +#' fit <- mkinfit(SFO_SFO, FOCUS_2006_D, quiet = TRUE) #' fit <- mkinfit(SFO_SFO, FOCUS_2006_D, quiet = TRUE, error_model = "tc") #' plot(fit) #' plot_res(fit) diff --git a/man/mkinfit.Rd b/man/mkinfit.Rd index b728aded..db3f5f3e 100644 --- a/man/mkinfit.Rd +++ b/man/mkinfit.Rd @@ -26,7 +26,6 @@ mkinfit( quiet = FALSE, atol = 1e-08, rtol = 1e-10, - n.outtimes = 10, error_model = c("const", "obs", "tc"), error_model_algorithm = c("auto", "d_3", "direct", "twostep", "threestep", "fourstep", "IRLS", "OLS"), @@ -141,11 +140,6 @@ is 1e-8, lower than in \code{\link{lsoda}}.} \item{rtol}{Absolute error tolerance, passed to \code{\link{ode}}. Default is 1e-10, much lower than in \code{\link{lsoda}}.} -\item{n.outtimes}{The length of the dataseries that is produced by the model -prediction function \code{\link{mkinpredict}}. This impacts the accuracy -of the numerical solver if that is used (see \code{solution_type} -argument.} - \item{error_model}{If the error model is "const", a constant standard deviation is assumed. diff --git a/man/mkinmod.Rd b/man/mkinmod.Rd index 2ba917d6..a5736be7 100644 --- a/man/mkinmod.Rd +++ b/man/mkinmod.Rd @@ -82,7 +82,7 @@ The IORE submodel is not well tested for metabolites. When using this \examples{ # Specify the SFO model (this is not needed any more, as we can now mkinfit("SFO", ...) -SFO <- mkinmod(parent = list(type = "SFO")) +SFO <- mkinmod(parent = mkinsub("SFO")) # One parent compound, one metabolite, both single first order SFO_SFO <- mkinmod( diff --git a/man/mkinpredict.Rd b/man/mkinpredict.Rd index f7e4acfc..1ba5b5f8 100644 --- a/man/mkinpredict.Rd +++ b/man/mkinpredict.Rd @@ -91,7 +91,7 @@ FALSE).} solver is used.} } \value{ -A matrix in the same format as the output of \code{\link{ode}}. +A data frame with the solution in wide format } \description{ This function produces a time series for all the observed variables in a @@ -131,30 +131,35 @@ mkinpredict(SFO, c(k_degradinol = 0.3), c(degradinol = 100), # Check compiled model versions - they are faster than the eigenvalue based solutions! SFO_SFO = mkinmod(parent = list(type = "SFO", to = "m1"), - m1 = list(type = "SFO"), use_of_ff = "min") + m1 = list(type = "SFO"), use_of_ff = "max") if(require(rbenchmark)) { - benchmark( - eigen = mkinpredict(SFO_SFO, c(k_parent_m1 = 0.05, k_parent_sink = 0.1, k_m1_sink = 0.01), + benchmark(replications = 10, order = "relative", columns = c("test", "relative", "elapsed"), + eigen = mkinpredict(SFO_SFO, + c(k_parent = 0.15, f_parent_to_m1 = 0.5, k_m1 = 0.01), c(parent = 100, m1 = 0), seq(0, 20, by = 0.1), solution_type = "eigen")[201,], - deSolve_compiled = mkinpredict(SFO_SFO, - c(k_parent_m1 = 0.05, k_parent_sink = 0.1, k_m1_sink = 0.01), + deSolve_compiled = mkinpredict(SFO_SFO, + c(k_parent = 0.15, f_parent_to_m1 = 0.5, k_m1 = 0.01), c(parent = 100, m1 = 0), seq(0, 20, by = 0.1), solution_type = "deSolve")[201,], - deSolve = mkinpredict(SFO_SFO, c(k_parent_m1 = 0.05, k_parent_sink = 0.1, k_m1_sink = 0.01), - c(parent = 100, m1 = 0), seq(0, 20, by = 0.1), - solution_type = "deSolve", use_compiled = FALSE)[201,], - replications = 10) + deSolve = mkinpredict(SFO_SFO, + c(k_parent = 0.15, f_parent_to_m1 = 0.5, k_m1 = 0.01), + c(parent = 100, m1 = 0), seq(0, 20, by = 0.1), + solution_type = "deSolve", use_compiled = FALSE)[201,], + analytical = mkinpredict(SFO_SFO, + c(k_parent = 0.15, f_parent_to_m1 = 0.5, k_m1 = 0.01), + c(parent = 100, m1 = 0), seq(0, 20, by = 0.1), + solution_type = "analytical", use_compiled = FALSE)[201,]) } - -# Since mkin 0.9.49.11 we also have analytical solutions for some models, including SFO-SFO -# deSolve = mkinpredict(SFO_SFO, c(k_parent_m1 = 0.05, k_parent_sink = 0.1, k_m1_sink = 0.01), -# c(parent = 100, m1 = 0), seq(0, 20, by = 0.1), -# solution_type = "analytical", use_compiled = FALSE)[201,], + analytical = mkinpredict(SFO_SFO, + c(k_parent = 0.15, f_parent_to_m1 = 0.5, k_m1 = 0.01), + c(parent = 100, m1 = 0), seq(0, 20, by = 0.1), + solution_type = "analytical", use_compiled = FALSE)[201,] \dontrun{ # Predict from a fitted model f <- mkinfit(SFO_SFO, FOCUS_2006_C, quiet = TRUE) + f <- mkinfit(SFO_SFO, FOCUS_2006_C, quiet = TRUE, solution_type = "analytical") head(mkinpredict(f)) } diff --git a/man/plot.mkinfit.Rd b/man/plot.mkinfit.Rd index c3f3134a..4235557e 100644 --- a/man/plot.mkinfit.Rd +++ b/man/plot.mkinfit.Rd @@ -144,6 +144,7 @@ latex is being used for the formatting of the chi2 error level, if \dontrun{ SFO_SFO <- mkinmod(parent = mkinsub("SFO", "m1", full = "Parent"), m1 = mkinsub("SFO", full = "Metabolite M1" )) +fit <- mkinfit(SFO_SFO, FOCUS_2006_D, quiet = TRUE) fit <- mkinfit(SFO_SFO, FOCUS_2006_D, quiet = TRUE, error_model = "tc") plot(fit) plot_res(fit) diff --git a/test.log b/test.log index 489550a2..7788f109 100644 --- a/test.log +++ b/test.log @@ -2,34 +2,35 @@ Loading mkin Testing mkin ✔ | OK F W S | Context ✔ | 2 | Export dataset for reading into CAKE -✔ | 13 | Results for FOCUS D established in expertise for UBA (Ranke 2014) [3.1 s] -✔ | 4 | Calculation of FOCUS chi2 error levels [1.9 s] -✔ | 7 | Fitting the SFORB model [9.6 s] +✔ | 13 | Results for FOCUS D established in expertise for UBA (Ranke 2014) [3.3 s] +✔ | 4 | Calculation of FOCUS chi2 error levels [2.0 s] +✔ | 7 | Fitting the SFORB model [9.1 s] +✔ | 1 | Analytical solutions for coupled models [1.1 s] ✔ | 5 | Calculation of Akaike weights -✔ | 10 | Confidence intervals and p-values [8.3 s] -✔ | 14 | Error model fitting [22.1 s] -✔ | 6 | Test fitting the decline of metabolites from their maximum [0.7 s] -✔ | 1 | Fitting the logistic model [0.7 s] +✔ | 10 | Confidence intervals and p-values [1.0 s] +✔ | 14 | Error model fitting [3.7 s] +✔ | 6 | Test fitting the decline of metabolites from their maximum [0.2 s] +✔ | 1 | Fitting the logistic model [0.2 s] ✔ | 1 | Test dataset class mkinds used in gmkin -✔ | 12 | Special cases of mkinfit calls [2.1 s] +✔ | 12 | Special cases of mkinfit calls [0.6 s] ✔ | 8 | mkinmod model generation and printing [0.2 s] ✔ | 3 | Model predictions with mkinpredict [0.4 s] -✔ | 16 | Evaluations according to 2015 NAFTA guidance [3.8 s] -✔ | 9 | Nonlinear mixed-effects models [11.9 s] +✔ | 16 | Evaluations according to 2015 NAFTA guidance [1.4 s] +✔ | 9 | Nonlinear mixed-effects models [12.1 s] ✔ | 4 | Calculation of maximum time weighted average concentrations (TWAs) [2.4 s] ✔ | 3 | Summary -✔ | 14 | Plotting [4.0 s] +✔ | 14 | Plotting [3.5 s] ✔ | 4 | AIC calculation ✔ | 4 | Residuals extracted from mkinfit models -✔ | 2 | Complex test case from Schaefer et al. (2007) Piacenza paper [3.8 s] +✔ | 2 | Complex test case from Schaefer et al. (2007) Piacenza paper [4.2 s] ✔ | 1 | Summaries of old mkinfit objects -✔ | 4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [6.0 s] -✔ | 9 | Hypothesis tests [21.3 s] +✔ | 4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [6.4 s] +✔ | 9 | Hypothesis tests [17.8 s] ══ Results ═════════════════════════════════════════════════════════════════════ -Duration: 102.5 s +Duration: 69.9 s -OK: 156 +OK: 157 Failed: 0 Warnings: 0 Skipped: 0 diff --git a/tests/testthat/FOCUS_2006_D.csf b/tests/testthat/FOCUS_2006_D.csf index f27df430..f113a668 100644 --- a/tests/testthat/FOCUS_2006_D.csf +++ b/tests/testthat/FOCUS_2006_D.csf @@ -5,7 +5,7 @@ Description: MeasurementUnits: % AR TimeUnits: days Comments: Created using mkin::CAKE_export -Date: 2020-05-08 +Date: 2020-05-09 Optimiser: IRLS [Data] diff --git a/tests/testthat/test_analytical.R b/tests/testthat/test_analytical.R new file mode 100644 index 00000000..3d30e042 --- /dev/null +++ b/tests/testthat/test_analytical.R @@ -0,0 +1,11 @@ +context("Analytical solutions for coupled models") + +test_that("The analytical solution of SFO-SFO is correct", { + f_sfo_sfo.ff.analytical <- mkinfit(SFO_SFO.ff, + subset(FOCUS_2006_D, value != 0), + quiet = TRUE) + expect_equal( + parms(f_sfo_sfo.ff), + parms(f_sfo_sfo.ff.analytical) + ) +}) diff --git a/tests/testthat/test_tests.R b/tests/testthat/test_tests.R index 3d1c136c..b7c3c548 100644 --- a/tests/testthat/test_tests.R +++ b/tests/testthat/test_tests.R @@ -49,7 +49,7 @@ test_that("Updating fitted models works", { error_model = "tc", quiet = TRUE) f_soil_1_nw <- update(f_soil_1_tc, error_model = "const") - f_soil_1_nw_A2 <- update(f_soil_1_nw, fixed_parms = c(k_A2 = 0)) + f_soil_1_nw_A2 <- suppressWarnings(update(f_soil_1_nw, fixed_parms = c(k_A2 = 0))) test_nw_tc <- lrtest(f_soil_1_nw, f_soil_1_tc) expect_equivalent(test_nw_tc[["2", "Pr(>Chisq)"]], 2.113e-6) test_nw_A2 <- lrtest(f_soil_1_nw, f_soil_1_nw_A2) -- cgit v1.2.1