## ----------------------------------------------------------------------------- ## Ordinary Differential Equations Solver function. ## ----------------------------------------------------------------------------- # Some of the CAKE R modules are based on mkin. # # Modifications developed by Hybrid Intelligence (formerly Tessella), part of # Capgemini Engineering, for Syngenta, Copyright (C) 2011-2022 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 # 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)) }