From ae64167d93bfae36158578f0a1c7771e6bab9350 Mon Sep 17 00:00:00 2001
From: Johannes Ranke <jranke@uni-bremen.de>
Date: Thu, 4 Jun 2020 13:27:44 +0200
Subject: Version 3.4 as just publicly announced

Peter Rainbird just announced the release on the PFMODELS email list.
---
 CakeOdeSolve.R | 105 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 105 insertions(+)
 create mode 100644 CakeOdeSolve.R

(limited to 'CakeOdeSolve.R')

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
-- 
cgit v1.2.1