diff options
Diffstat (limited to 'CakeOdeSolve.R')
| -rw-r--r-- | CakeOdeSolve.R | 105 |
1 files changed, 105 insertions, 0 deletions
diff --git a/CakeOdeSolve.R b/CakeOdeSolve.R new file mode 100644 index 0000000..263e57c --- /dev/null +++ b/CakeOdeSolve.R @@ -0,0 +1,105 @@ +## ----------------------------------------------------------------------------- +## Ordinary Differential Equations Solver function. +## ----------------------------------------------------------------------------- +# Some of the CAKE R modules are based on mkin. +# +# Modifications developed by Tessella for Syngenta, Copyright (C) 2011-2020 Syngenta +# Tessella Project Reference: 6245, 7247, 8361, 7414, 10091 +# +# The CAKE R modules are free software: you can +# redistribute them and/or modify them under the +# terms of the GNU General Public License as published by the Free Software +# Foundation, either version 3 of the License, or (at your option) any later +# version. +# +# This program is distributed in the hope that it will be useful, but WITHOUT +# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +# details. +# +# You should have received a copy of the GNU General Public License along with +# this program. If not, see <http://www.gnu.org/licenses/> + + +# fit: Fit with initial (state) values and parameters for the ODE system. +# outtimes: Time sequence for output. +# solution: Whether to use analytical, eigenvectors, or general ode solver to solve the ODE system. +# atol: The tolerance to apply to the ODE solver. +CakeOdeSolve <- function(fit, outtimes, solution, atol) +{ + fixed <- fit$fixed$value + names(fixed) <- rownames(fit$fixed) + parms.all <- c(fit$par, fixed) + ininames <- c( + rownames(subset(fit$start, type == "state")), + rownames(subset(fit$fixed, type == "state"))) + odeini <- parms.all[ininames] + names(odeini) <- gsub("_0$", "", names(odeini)) + odenames <- c( + rownames(subset(fit$start, type == "deparm")), + rownames(subset(fit$fixed, type == "deparm"))) + odeparms <- parms.all[odenames] + odeini <- AdjustOdeInitialValues(odeini, fit, odeparms) + + evalparse <- function(string) + { + eval(parse(text = string), as.list(c(odeparms, odeini))) + } + + odeResult <- numeric() + + if (solution == "analytical") + { + parent.type = names(fit$map[[1]])[1] + parent.name = names(fit$diffs)[[1]] + ode <- switch(parent.type, + SFO = SFO.solution(outtimes, + evalparse(parent.name), + evalparse(paste("k", parent.name, sep = "_"))), + FOMC = FOMC.solution(outtimes, + evalparse(parent.name), + evalparse("alpha"), evalparse("beta")), + DFOP = DFOP.solution(outtimes, + evalparse(parent.name), + evalparse(paste("k1", parent.name, sep = "_")), + evalparse(paste("k2", parent.name, sep = "_")), + evalparse(paste("g", parent.name, sep = "_"))), + HS = HS.solution(outtimes, + evalparse(parent.name), + evalparse("k1"), evalparse("k2"), + evalparse("tb")), + IORE = IORE.solution(outtimes, + evalparse(parent.name), + evalparse(paste("k", parent.name, sep = "_")), + evalparse("N")) + ) + odeResult <- cbind(outtimes, ode) + dimnames(odeResult) <- list(outtimes, c("time", parent.name)) + } + else if (solution == "eigen") + { + coefmat.num <- matrix(sapply(as.vector(fit$coefmat), evalparse), + nrow = length(odeini)) + e <- eigen(coefmat.num) + c <- solve(e$vectors, odeini) + f.out <- function(t) { + e$vectors %*% diag(exp(e$values * t), nrow = length(odeini)) %*% c + } + ode <- matrix(mapply(f.out, outtimes), + nrow = length(odeini), ncol = length(outtimes)) + dimnames(ode) <- list(names(odeini), NULL) + odeResult <- cbind(time = outtimes, t(ode)) + } + else if (solution == "deSolve") + { + odeResult <- ode( + y = odeini, + times = outtimes, + func = fit$mkindiff, + parms = odeparms, + atol = atol + ) + } + + return(data.frame(odeResult)) +}
\ No newline at end of file |