1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
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))
}
|
