From ec574cff822a1238138c0aa69b3d1459bdc3dfa8 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Fri, 19 Jun 2015 17:46:11 +0200 Subject: Use odeintr instead of ccSolve for compiling models --- R/IORE.solution.R | 4 +-- R/endpoints.R | 2 +- R/mkinerrmin.R | 6 ++-- R/mkinfit.R | 23 ++++++------ R/mkinmod.R | 95 ++++++++++++++++++++++++++++++++------------------ R/mkinpredict.R | 50 +++++++++++++++----------- R/plot.mkinfit.R | 2 +- R/transform_odeparms.R | 12 +++---- 8 files changed, 115 insertions(+), 79 deletions(-) (limited to 'R') diff --git a/R/IORE.solution.R b/R/IORE.solution.R index 9546ce56..5405be96 100644 --- a/R/IORE.solution.R +++ b/R/IORE.solution.R @@ -1,4 +1,4 @@ -IORE.solution <- function(t, parent.0, k.iore, N) +IORE.solution <- function(t, parent.0, k__iore, N) { - parent = (parent.0^(1 - N) - (1 - N) * k.iore * t)^(1/(1 - N)) + parent = (parent.0^(1 - N) - (1 - N) * k__iore * t)^(1/(1 - N)) } diff --git a/R/endpoints.R b/R/endpoints.R index ae9aa3ec..001595bc 100644 --- a/R/endpoints.R +++ b/R/endpoints.R @@ -50,7 +50,7 @@ endpoints <- function(fit) { ep$distimes[obs_var, c("DT50back")] = DT50_back } if (type == "IORE") { - k_names = grep(paste("k.iore", obs_var, sep="_"), names(parms.all), value=TRUE) + k_names = grep(paste("k__iore", obs_var, sep="_"), names(parms.all), value=TRUE) k_tot = sum(parms.all[k_names]) # From the NAFTA kinetics guidance, p. 5 n = parms.all[paste("N", obs_var, sep = "_")] diff --git a/R/mkinerrmin.R b/R/mkinerrmin.R index b28235a6..6103dd1d 100644 --- a/R/mkinerrmin.R +++ b/R/mkinerrmin.R @@ -61,8 +61,8 @@ mkinerrmin <- function(fit, alpha = 0.05) n.k.optim <- length(grep(paste("^k", obs_var, sep="_"), names(parms.optim))) n.k.optim <- n.k.optim + length(grep(paste("^log_k", obs_var, sep="_"), names(parms.optim))) - n.k.iore.optim <- length(grep(paste("^k.iore", obs_var, sep="_"), names(parms.optim))) - n.k.iore.optim <- n.k.iore.optim + length(grep(paste("^log_k.iore", obs_var, + n.k__iore.optim <- length(grep(paste("^k__iore", obs_var, sep="_"), names(parms.optim))) + n.k__iore.optim <- n.k__iore.optim + length(grep(paste("^log_k__iore", obs_var, sep = "_"), names(parms.optim))) @@ -82,7 +82,7 @@ mkinerrmin <- function(fit, alpha = 0.05) } } - n.optim <- sum(n.initials.optim, n.k.optim, n.k.iore.optim, n.N.optim, n.ff.optim) + n.optim <- sum(n.initials.optim, n.k.optim, n.k__iore.optim, n.N.optim, n.ff.optim) # FOMC, DFOP and HS parameters are only counted if we are looking at the # first variable in the model which is always the source variable diff --git a/R/mkinfit.R b/R/mkinfit.R index 5118519a..ce380243 100644 --- a/R/mkinfit.R +++ b/R/mkinfit.R @@ -26,9 +26,8 @@ mkinfit <- function(mkinmod, observed, state.ini = "auto", fixed_parms = NULL, fixed_initials = names(mkinmod$diffs)[-1], - solution_type = c("auto", "analytical", "eigen", "deSolve"), + solution_type = c("auto", "analytical", "eigen", "deSolve", "odeintr"), method.ode = "lsoda", - use_compiled = "auto", method.modFit = c("Port", "Marq", "SANN", "Nelder-Mead", "BFGS", "CG", "L-BFGS-B"), maxit.modFit = "auto", control.modFit = list(), @@ -208,8 +207,7 @@ mkinfit <- function(mkinmod, observed, # Decide if the solution of the model can be based on a simple analytical # formula, the spectral decomposition of the matrix (fundamental system) # or a numeric ode solver from the deSolve package - # Prefer deSolve over eigen if a compiled model is present and use_compiled - # is not set to FALSE + # Prefer odeintr if a compiled model is present solution_type = match.arg(solution_type) if (solution_type == "analytical" && length(mkinmod$spec) > 1) stop("Analytical solution not implemented for models with metabolites.") @@ -219,8 +217,11 @@ mkinfit <- function(mkinmod, observed, if (length(mkinmod$spec) == 1) { solution_type = "analytical" } else { - if (!is.null(mkinmod$compiled) & use_compiled[1] != FALSE) { - solution_type = "deSolve" + if (!is.null(mkinmod$e$m)) { + solution_type = "odeintr" + if (!missing(atol) | !missing(rtol)) { + message("When fitting with solution type odeintr, atol and rtol can only be set in the call to mkinmod()") + } } else { if (is.matrix(mkinmod$coefmat)) { solution_type = "eigen" @@ -269,9 +270,8 @@ mkinfit <- function(mkinmod, observed, out <- mkinpredict(mkinmod, parms, odeini, outtimes, solution_type = solution_type, - use_compiled = use_compiled, method.ode = method.ode, - atol = atol, rtol = rtol, ...) + atol.deSolve = atol, rtol.deSolve = rtol, ...) assign("out_predicted", out, inherits=TRUE) @@ -289,9 +289,8 @@ mkinfit <- function(mkinmod, observed, out_plot <- mkinpredict(mkinmod, parms, odeini, outtimes_plot, solution_type = solution_type, - use_compiled = use_compiled, method.ode = method.ode, - atol = atol, rtol = rtol, ...) + atol.deSolve = atol, rtol.deSolve = rtol, ...) plot(0, type="n", xlim = range(observed$time), ylim = c(0, max(observed$value, na.rm=TRUE)), @@ -319,8 +318,8 @@ mkinfit <- function(mkinmod, observed, if (!transform_rates) { index_k <- grep("^k_", names(lower)) lower[index_k] <- 0 - index_k.iore <- grep("^k.iore_", names(lower)) - lower[index_k.iore] <- 0 + index_k__iore <- grep("^k__iore_", names(lower)) + lower[index_k__iore] <- 0 other_rate_parms <- intersect(c("alpha", "beta", "k1", "k2", "tb"), names(lower)) lower[other_rate_parms] <- 0 } diff --git a/R/mkinmod.R b/R/mkinmod.R index e865fa27..4991be22 100644 --- a/R/mkinmod.R +++ b/R/mkinmod.R @@ -16,12 +16,17 @@ # You should have received a copy of the GNU General Public License along with # this program. If not, see }}} -mkinmod <- function(..., use_of_ff = "min", speclist = NULL, quiet = FALSE) +mkinmod <- function(..., use_of_ff = "min", speclist = NULL, quiet = FALSE, + odeintr_compile = c("auto", "no", "yes"), + odeintr_method = "rk5_i", + odeintr_atol = 1e-8, odeintr_rtol = 1e-10) { if (is.null(speclist)) spec <- list(...) else spec <- speclist obs_vars <- names(spec) + odeintr_compile <- match.arg(odeintr_compile) + # Check if any of the names of the observed variables contains any other for (obs_var in obs_vars) { if (length(grep(obs_var, obs_vars)) > 1) stop("Sorry, variable names can not contain each other") @@ -99,7 +104,7 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL, quiet = FALSE) # If sink is required, add first-order/IORE sink term k_compound_sink <- paste("k", box_1, "sink", sep = "_") if(spec[[varname]]$type == "IORE") { - k_compound_sink <- paste("k.iore", box_1, "sink", sep = "_") + k_compound_sink <- paste("k__iore", box_1, "sink", sep = "_") } parms <- c(parms, k_compound_sink) decline_term <- paste(k_compound_sink, "*", box_1) @@ -114,7 +119,7 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL, quiet = FALSE) } else { k_compound <- paste("k", box_1, sep = "_") if(spec[[varname]]$type == "IORE") { - k_compound <- paste("k.iore", box_1, sep = "_") + k_compound <- paste("k__iore", box_1, sep = "_") } parms <- c(parms, k_compound) decline_term <- paste(k_compound, "*", box_1) @@ -127,17 +132,18 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL, quiet = FALSE) } #}}} if(spec[[varname]]$type == "FOMC") { # {{{ Add FOMC decline term # From p. 53 of the FOCUS kinetics report, without the power function so it works in C - decline_term <- paste("(alpha/beta) * 1/((time/beta) + 1) *", box_1) + decline_term <- paste("(alpha/beta) * 1/(( time /beta) + 1) *", box_1) parms <- c(parms, "alpha", "beta") } #}}} if(spec[[varname]]$type == "DFOP") { # {{{ Add DFOP decline term # From p. 57 of the FOCUS kinetics report - decline_term <- paste("((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 * time)) / (g * exp(-k1 * time) + (1 - g) * exp(-k2 * time))) *", box_1) + decline_term <- paste("((k1 * g * exp(-k1 * time ) + k2 * (1 - g) * exp(-k2 * time )) / (g * exp(-k1 * time ) + (1 - g) * exp(-k2 * time ))) *", box_1) parms <- c(parms, "k1", "k2", "g") } #}}} + HS_decline <- "ifelse(time <= tb, k1, k2)" # Used below for automatic translation to C++ if(spec[[varname]]$type == "HS") { # {{{ Add HS decline term # From p. 55 of the FOCUS kinetics report - decline_term <- paste("ifelse(time <= tb, k1, k2)", "*", box_1) + decline_term <- paste(HS_decline, "*", box_1) parms <- c(parms, "k1", "k2", "tb") } #}}} # Add origin decline term to box 1 (usually the only box, unless type is SFORB)#{{{ @@ -272,35 +278,58 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL, quiet = FALSE) model$coefmat <- m }#}}} - # Create a function compiled from C code if more than one observed variable and ccSolve is available #{{{ - if (length(obs_vars) > 1 & requireNamespace("ccSolve", quietly = TRUE)) { - diffs.C <- paste(diffs, collapse = ";\n") - diffs.C <- paste0(diffs.C, ";") - for (i in seq_along(diffs)) { - obs_var <- names(diffs)[i] + # Try to create a function compiled from C code if more than one observed variable and g++ can be called #{{{ + if (length(obs_vars) > 1) { + try(gcc_version <- system("g++ -v 2>&1", intern = TRUE)) + if (!inherits(gcc_version, "try-error") & odeintr_compile != "no") { - # Replace d_... terms by f[i-1] - # First line - pattern <- paste0("^d_", obs_var) - replacement <- paste0("\nf[", i - 1, "]") - diffs.C <- gsub(pattern, replacement, diffs.C) - # Other lines - pattern <- paste0("\\nd_", obs_var) - replacement <- paste0("\nf[", i - 1, "]") - diffs.C <- gsub(pattern, replacement, diffs.C) + if ((odeintr_compile == "auto" & mat == FALSE) | odeintr_compile == "yes") { + # HS + diffs <- gsub(HS_decline, "(t <= tb ? k1 : k2)", diffs, fixed = TRUE) - # Replace names of observed variables by y[i], - # making the implicit assumption that the observed variables only occur after "* " - pattern <- paste0("\\* ", obs_var) - replacement <- paste0("* y[", i - 1, "]") - diffs.C <- gsub(pattern, replacement, diffs.C) - } - if (sum(sapply(spec, function(x) x$type %in% - c("SFO", "FOMC", "DFOP", "SFORB"))) == length(spec)) { - if (!quiet) message("Compiling differential equation model from auto-generated C code...") - model$compiled <- ccSolve::compile.ode(diffs.C, language = "C", - parms = parms, - declaration = "double time = *t;") + diffs.C <- paste(diffs, collapse = ";\n") + diffs.C <- paste0(diffs.C, ";") + + for (i in seq_along(diffs)) { + state_var <- names(diffs)[i] + + # IORE + if (state_var %in% obs_vars) { + if (spec[[state_var]]$type == "IORE") { + diffs.C <- gsub(paste0(state_var, "^N_", state_var), + paste0("pow(x[", i - 1, "], N_", state_var, ")"), + diffs.C, fixed = TRUE) + } + } + + # Replace d_... terms by dxdt[i-1] + # First line + pattern <- paste0("^d_", state_var) + replacement <- paste0("\ndxdt[", i - 1, "]") + diffs.C <- gsub(pattern, replacement, diffs.C) + # Other lines + pattern <- paste0("\\nd_", state_var) + replacement <- paste0("\ndxdt[", i - 1, "]") + diffs.C <- gsub(pattern, replacement, diffs.C) + + # Replace names of observed variables by x[i], + # making the implicit assumption that the observed variables only occur after "* " + pattern <- paste0("\\* ", state_var) + replacement <- paste0("* x[", i - 1, "]") + diffs.C <- gsub(pattern, replacement, diffs.C) + + # Handle time, which needs to be enclosed by spaces (except for HS which is handled separately) + diffs.C <- gsub(" time ", " t ", diffs.C) + } + if (!quiet) message("Compiling differential equation model from auto-generated C++ code...") + model$e <- new.env() + model$e$code <- try(compile_sys("m", diffs.C, pars = parms, method = odeintr_method, + atol = odeintr_atol, rtol = odeintr_rtol, env = model$e)) + + if (inherits(model$e$code, "try-error")) { + if (!quiet) message("Compilation failed, model solution will be slower than necessary") + } + } } } # }}} diff --git a/R/mkinpredict.R b/R/mkinpredict.R index 3121f1d7..d91e87df 100644 --- a/R/mkinpredict.R +++ b/R/mkinpredict.R @@ -18,11 +18,12 @@ # this program. If not, see mkinpredict <- function(mkinmod, odeparms, odeini, - outtimes, solution_type = "deSolve", - use_compiled = "auto", - method.ode = "lsoda", atol = 1e-8, rtol = 1e-10, + outtimes, solution_type = c("deSolve", "analytical", "eigen", "odeintr"), + method.ode = "lsoda", atol.deSolve = 1e-8, rtol.deSolve = 1e-10, map_output = TRUE, ...) { + solution_type = match.arg(solution_type) + # Get the names of the state variables in the model mod_vars <- names(mkinmod$diffs) @@ -54,8 +55,8 @@ mkinpredict <- function(mkinmod, odeparms, odeini, IORE = IORE.solution(outtimes, evalparse(parent.name), ifelse(mkinmod$use_of_ff == "min", - evalparse(paste("k.iore", parent.name, "sink", sep="_")), - evalparse(paste("k.iore", parent.name, sep="_"))), + evalparse(paste("k__iore", parent.name, "sink", sep="_")), + evalparse(paste("k__iore", parent.name, sep="_"))), evalparse("N_parent")), DFOP = DFOP.solution(outtimes, evalparse(parent.name), @@ -88,30 +89,26 @@ mkinpredict <- function(mkinmod, odeparms, odeini, names(out) <- c("time", mod_vars) } if (solution_type == "deSolve") { - if (!is.null(mkinmod$compiled) & use_compiled[1] != FALSE) { - mkindiff <- mkinmod$compiled - } else { - mkindiff <- function(t, state, parms) { + mkindiff <- function(t, state, parms) { - time <- t - diffs <- vector() - for (box in names(mkinmod$diffs)) - { - diffname <- paste("d", box, sep="_") - diffs[diffname] <- with(as.list(c(time, state, parms)), - eval(parse(text=mkinmod$diffs[[box]]))) - } - return(list(c(diffs))) + time <- t + diffs <- vector() + for (box in names(mkinmod$diffs)) + { + diffname <- paste("d", box, sep="_") + diffs[diffname] <- with(as.list(c(time, state, parms)), + eval(parse(text=mkinmod$diffs[[box]]))) } - } + return(list(c(diffs))) + } out <- ode( y = odeini, times = outtimes, func = mkindiff, parms = odeparms[mkinmod$parms], # Order matters when using compiled models method = method.ode, - atol = atol, - rtol = rtol, + atol = atol.deSolve, + rtol = rtol.deSolve, ... ) if (sum(is.na(out)) > 0) { @@ -119,6 +116,17 @@ mkinpredict <- function(mkinmod, odeparms, odeini, "NaN values occurred in output from ode()") } } + if (solution_type == "odeintr") { + if (is.null(mkinmod$e$m)) stop("Method odeintr can not be used as no compiled version of the model is available") + odeparms_argstring = "" + for (parname in names(odeparms)) { + odeparms_argstring = paste0(odeparms_argstring, parname, " = ", odeparms[parname], ", ") + } + odeparms_argstring = gsub(", $", "", odeparms_argstring) + with(as.list(odeparms_argstring), eval(parse(text=paste0("mkinmod$e$m_set_params(", odeparms_argstring, ")")))) + out <- mkinmod$e$m_at(init = odeini, times = outtimes) + names(out) <- c("time", names(mkinmod$diffs)) + } if (map_output) { # Output transformation for models with unobserved compartments like SFORB out_mapped <- data.frame(time = out[,"time"]) diff --git a/R/plot.mkinfit.R b/R/plot.mkinfit.R index 31746fb8..1523be7d 100644 --- a/R/plot.mkinfit.R +++ b/R/plot.mkinfit.R @@ -55,7 +55,7 @@ plot.mkinfit <- function(x, fit = x, odeparms <- parms.all[odenames] out <- mkinpredict(fit$mkinmod, odeparms, odeini, outtimes, - solution_type = solution_type, atol = fit$atol, rtol = fit$rtol) + solution_type = solution_type, atol.deSolve = fit$atol, rtol.deSolve = fit$rtol) # Set up the plot if not to be added to an existing plot if (add == FALSE) { diff --git a/R/transform_odeparms.R b/R/transform_odeparms.R index aa70a0b3..0946ff14 100644 --- a/R/transform_odeparms.R +++ b/R/transform_odeparms.R @@ -33,8 +33,8 @@ transform_odeparms <- function(parms, mkinmod, # Log transformation for rate constants if requested k <- parms[grep("^k_", names(parms))] - k.iore <- parms[grep("^k.iore_", names(parms))] - k <- c(k, k.iore) + k__iore <- parms[grep("^k__iore_", names(parms))] + k <- c(k, k__iore) if (length(k) > 0) { if(transform_rates) { transparms[paste0("log_", names(k))] <- log(k) @@ -113,8 +113,8 @@ backtransform_odeparms <- function(transparms, mkinmod, # Exponential transformation for rate constants if(transform_rates) { trans_k <- transparms[grep("^log_k_", names(transparms))] - trans_k.iore <- transparms[grep("^log_k.iore_", names(transparms))] - trans_k = c(trans_k, trans_k.iore) + trans_k__iore <- transparms[grep("^log_k__iore_", names(transparms))] + trans_k = c(trans_k, trans_k__iore) if (length(trans_k) > 0) { k_names <- gsub("^log_k", "k", names(trans_k)) parms[k_names] <- exp(trans_k) @@ -122,8 +122,8 @@ backtransform_odeparms <- function(transparms, mkinmod, } else { trans_k <- transparms[grep("^k_", names(transparms))] parms[names(trans_k)] <- trans_k - trans_k.iore <- transparms[grep("^k.iore_", names(transparms))] - parms[names(trans_k.iore)] <- trans_k.iore + trans_k__iore <- transparms[grep("^k__iore_", names(transparms))] + parms[names(trans_k__iore)] <- trans_k__iore } # Do not transform exponents in IORE models -- cgit v1.2.1