summaryrefslogtreecommitdiff
path: root/CakeOdeSolve.R
diff options
context:
space:
mode:
authorJohannes Ranke <jranke@uni-bremen.de>2020-06-04 13:27:44 +0200
committerJohannes Ranke <jranke@uni-bremen.de>2020-06-04 13:27:44 +0200
commitae64167d93bfae36158578f0a1c7771e6bab9350 (patch)
treea29ed6dc384956d6b35587c628f8ff035e09c327 /CakeOdeSolve.R
parent1684a82b15dee35812c1340e26d721ee64a28751 (diff)
Version 3.4 as just publicly announcedv3.4
Peter Rainbird just announced the release on the PFMODELS email list.
Diffstat (limited to 'CakeOdeSolve.R')
-rw-r--r--CakeOdeSolve.R105
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

Contact - Imprint