diff options
Diffstat (limited to 'R')
| -rw-r--r-- | R/mkinfit.R | 6 | ||||
| -rw-r--r-- | R/mkinmod.R | 40 | ||||
| -rw-r--r-- | R/mkinpredict.R | 152 |
3 files changed, 104 insertions, 94 deletions
diff --git a/R/mkinfit.R b/R/mkinfit.R index 5c092612..2b7e71cb 100644 --- a/R/mkinfit.R +++ b/R/mkinfit.R @@ -185,13 +185,13 @@ if(getRversion() >= '2.15.1') utils::globalVariables(c("name", "time", "value")) #' # Fit the model to the FOCUS example dataset D using defaults #' print(system.time(fit <- mkinfit(SFO_SFO, FOCUS_2006_D, #' solution_type = "eigen", quiet = TRUE))) -#' coef(fit) +#' parms(fit) #' endpoints(fit) #' \dontrun{ #' # deSolve is slower when no C compiler (gcc) was available during model generation #' print(system.time(fit.deSolve <- mkinfit(SFO_SFO, FOCUS_2006_D, #' solution_type = "deSolve"))) -#' coef(fit.deSolve) +#' parms(fit.deSolve) #' endpoints(fit.deSolve) #' } #' @@ -926,6 +926,6 @@ mkinfit <- function(mkinmod, observed, fit$version <- as.character(utils::packageVersion("mkin")) fit$Rversion <- paste(R.version$major, R.version$minor, sep=".") - class(fit) <- c("mkinfit", "modFit") + class(fit) <- c("mkinfit") return(fit) } diff --git a/R/mkinmod.R b/R/mkinmod.R index ca1402fd..099e1155 100644 --- a/R/mkinmod.R +++ b/R/mkinmod.R @@ -421,6 +421,46 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL, quiet = FALSE, verb } # }}} + # Attach a degradation function if an analytical solution is available {{{ + parent_type = spec[[1]]$type + parent_name = names(spec)[[1]] + if (length(spec) == 1) { + odeparm_map <- switch(parent_type, + SFO = c( + k = if(use_of_ff == "min") paste("k", parent_name, "sink", sep = "_") + else paste("k", parent_name, sep = "_")), + FOMC = c(alpha = "alpha", beta = "beta"), + IORE = c( + k__iore = if(use_of_ff == "min") paste("k__iore", parent_name, "sink", sep = "_") + else paste("k__iore", parent_name, sep = "_"), + N = paste("N", parent_name, sep = "_")), + DFOP = c(k1 = "k1", k2 = "k2", g = "g"), + HS = c(k1 = "k1", k2 = "k2", tb = "tb"), + SFORB = c( + k_12 = paste("k", parent_name, "free_bound", sep = "_"), + k_21 = paste("k", parent_name, "bound_free", sep = "_"), + k_1output = paste("k", parent_name, "free_sink", sep = "_")), + logistic = c(kmax = "kmax", k0 = "k0", r = "r") + ) + odeparm_rev_map <- names(odeparm_map) + names(odeparm_rev_map) <- odeparm_map + + model$deg_func <- function(odeini, odeparms, outtimes) { + parent_func <- getFromNamespace(paste0(parent_type, ".solution"), "mkin") + odeparm_list <- as.list(odeparms) + names(odeparm_list) <- odeparm_rev_map[names(odeparm_list)] + + values <- do.call(parent_func, + args = c( + list(t = outtimes, parent.0 = odeini[1]), + odeparm_list)) + out <- data.frame(outtimes, values) + names(out) <- c("time", parent_name) + return(out) + } + } + # }}} + class(model) <- "mkinmod" return(model) } diff --git a/R/mkinpredict.R b/R/mkinpredict.R index 16ee7903..0f8e83bb 100644 --- a/R/mkinpredict.R +++ b/R/mkinpredict.R @@ -41,57 +41,63 @@ #' @author Johannes Ranke #' @examples #' -#' SFO <- mkinmod(degradinol = mkinsub("SFO")) -#' # Compare solution types -#' mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), 0:20, -#' solution_type = "analytical") -#' mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), 0:20, -#' solution_type = "deSolve") -#' mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), 0:20, -#' solution_type = "deSolve", use_compiled = FALSE) -#' mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), 0:20, -#' solution_type = "eigen") +#' SFO <- mkinmod(degradinol = mkinsub("SFO")) +#' # Compare solution types +#' mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), 0:20, +#' solution_type = "analytical") +#' mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), 0:20, +#' solution_type = "deSolve") +#' mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), 0:20, +#' solution_type = "deSolve", use_compiled = FALSE) +#' mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), 0:20, +#' solution_type = "eigen") #' +#' # Compare integration methods to analytical solution +#' mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), 0:20, +#' solution_type = "analytical")[21,] +#' mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), 0:20, +#' method = "lsoda")[21,] +#' mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), 0:20, +#' method = "ode45")[21,] +#' mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), 0:20, +#' method = "rk4")[21,] +#' # rk4 is not as precise here #' -#' # Compare integration methods to analytical solution -#' mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), 0:20, -#' solution_type = "analytical")[21,] -#' mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), 0:20, -#' method = "lsoda")[21,] -#' mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), 0:20, -#' method = "ode45")[21,] -#' mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), 0:20, -#' method = "rk4")[21,] -#' # rk4 is not as precise here +#' # The number of output times used to make a lot of difference until the +#' # default for atol was adjusted +#' mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), +#' seq(0, 20, by = 0.1))[201,] +#' mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), +#' seq(0, 20, by = 0.01))[2001,] #' -#' # The number of output times used to make a lot of difference until the -#' # default for atol was adjusted -#' mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), -#' seq(0, 20, by = 0.1))[201,] -#' mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), -#' seq(0, 20, by = 0.01))[2001,] +#' # 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")) +#' if(require(rbenchmark)) { +#' benchmark( +#' eigen = 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 = "eigen")[201,], +#' deSolve_compiled = 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")[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) +#' } #' -#' # 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")) -#' system.time( -#' print(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 = "eigen")[201,])) -#' system.time( -#' print(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")[201,])) -#' system.time( -#' print(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,])) +#' # 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,], #' -#' \dontrun{ -#' # Predict from a fitted model -#' f <- mkinfit(SFO_SFO, FOCUS_2006_C) -#' head(mkinpredict(f)) -#' } +#' \dontrun{ +#' # Predict from a fitted model +#' f <- mkinfit(SFO_SFO, FOCUS_2006_C, quiet = TRUE) +#' head(mkinpredict(f)) +#' } #' #' @export mkinpredict <- function(x, odeparms, odeini, @@ -124,54 +130,17 @@ mkinpredict.mkinmod <- function(x, odeini <- odeini[mod_vars] } - # Create function for evaluation of expressions with ode parameters and initial values - evalparse <- function(string) - { - eval(parse(text=string), as.list(c(odeparms, odeini))) - } - - # Create a function calculating the differentials specified by the model - # if necessary if (solution_type == "analytical") { - parent.type = names(x$map[[1]])[1] - parent.name = names(x$diffs)[[1]] - o <- switch(parent.type, - SFO = SFO.solution(outtimes, - evalparse(parent.name), - ifelse(x$use_of_ff == "min", - evalparse(paste("k", parent.name, "sink", sep="_")), - evalparse(paste("k", parent.name, sep="_")))), - FOMC = FOMC.solution(outtimes, - evalparse(parent.name), - evalparse("alpha"), evalparse("beta")), - IORE = IORE.solution(outtimes, - evalparse(parent.name), - ifelse(x$use_of_ff == "min", - 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), - evalparse("k1"), evalparse("k2"), - evalparse("g")), - HS = HS.solution(outtimes, - evalparse(parent.name), - evalparse("k1"), evalparse("k2"), - evalparse("tb")), - SFORB = SFORB.solution(outtimes, - evalparse(parent.name), - evalparse(paste("k", parent.name, "bound", sep="_")), - evalparse(paste("k", sub("free", "bound", parent.name), "free", sep="_")), - evalparse(paste("k", parent.name, "sink", sep="_"))), - logistic = logistic.solution(outtimes, - evalparse(parent.name), - evalparse("kmax"), evalparse("k0"), - evalparse("r")) - ) - out <- data.frame(outtimes, o) - names(out) <- c("time", sub("_free", "", parent.name)) + out <- x$deg_func(odeini = odeini, + odeparms = odeparms, outtimes = outtimes) } + if (solution_type == "eigen") { + + evalparse <- function(string) { + eval(parse(text=string), as.list(c(odeparms, odeini))) + } + coefmat.num <- matrix(sapply(as.vector(x$coefmat), evalparse), nrow = length(mod_vars)) e <- eigen(coefmat.num) @@ -184,6 +153,7 @@ mkinpredict.mkinmod <- function(x, out <- data.frame(outtimes, t(o)) names(out) <- c("time", mod_vars) } + if (solution_type == "deSolve") { if (!is.null(x$cf) & use_compiled[1] != FALSE) { out <- ode( |