summaryrefslogblamecommitdiff
path: root/CakeOdeSolve.R
blob: f1c83f721fe2d9648c151aa9a80f07fd69b6c55f (plain) (tree) 
ae64167





fd9ea46


ae64167



































































































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
106





                                                                                

                                                                             


































































































                                                                                                   
## -----------------------------------------------------------------------------
## 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 <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))
}

Contact - Imprint