Version 1.4v1.4

1 files changed, 261 insertions, 0 deletions

diff --git a/CakeModel.R b/CakeModel.R
new file mode 100644
index 0000000..a93d7af
--- /dev/null
+++ b/CakeModel.R
@@ -0,0 +1,261 @@
+# Was Id: mkinmod.R 71 2010-09-12 01:13:36Z jranke 
+# $Id: CakeModel.R 216 2011-07-05 14:35:03Z nelr $
+
+# The CAKE R modules are based on mkin
+# Portions Johannes Ranke, 2010
+# Contact: mkin-devel@lists.berlios.de
+
+# This version has been modified to parameterise SFO as k and flow fractions
+# Modifications developed by Tessella Plc for Syngenta: Copyright (C) 2011  Syngenta
+# Tessella Project Reference: 6245
+
+#    This program is free software: you can redistribute it and/or modify
+#    it 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/>.”
+
+CakeModel <- function(...)
+{
+  spec <- list(...)
+  obs_vars <- names(spec)
+
+  # The returned model will be a list of character vectors, containing
+  # differential equations, parameter names and a mapping from model variables
+  # to observed variables. If possible, a matrix representation of the 
+  # differential equations is included
+  parms <- vector()
+  diffs <- vector()
+  map <- list()
+  if(spec[[1]]$type %in% c("FOMC", "DFOP", "HS")) {
+    mat = FALSE 
+    if(!is.null(spec[[1]]$to)) {
+      message <- paste("Only constant formation fractions over time are implemented.",
+        "Depending on the reason for the time dependence of degradation, this may be unrealistic",
+        sep="\n")
+      warning(message)
+    } else message <- "ok"
+  } else mat = TRUE
+  
+  mat = FALSE # XYZZY Not implemented yet assuming it should be
+  nlt = list()
+
+  # Establish list of differential equations
+  for (varname in obs_vars)
+  {
+    if(is.null(spec[[varname]]$type)) stop(
+      "Every argument to mkinmod must be a list containing a type component")
+    if(!spec[[varname]]$type %in% c("SFO", "FOMC", "DFOP", "HS", "SFORB")) stop(
+      "Available types are SFO, FOMC, DFOP, HS and SFORB only")
+    new_parms <- vector()
+
+    # New (sub)compartments (boxes) needed for the model type
+    new_boxes <- switch(spec[[varname]]$type,
+      SFO = varname,
+      FOMC = varname,
+      DFOP = varname,
+      HS = varname,
+      SFORB = paste(varname, c("free", "bound"), sep="_")
+    )
+    map[[varname]] <- new_boxes
+    names(map[[varname]]) <- rep(spec[[varname]]$type, length(new_boxes))
+
+    # Start a new differential equation for each new box
+    new_diffs <- paste("d_", new_boxes, " =", sep="")
+
+    # Turn on sink if not specified otherwise
+    if(is.null(spec[[varname]]$sink)) spec[[varname]]$sink <- TRUE
+    
+    #@@@ADD SFO k HERE !!!!!!!!!!!!!
+    # Construct and add SFO term and add SFO parameters if needed
+    if(spec[[varname]]$type == "SFO") {
+      # From p. 53 of the FOCUS kinetics report
+      k_term <- paste("k", new_boxes[[1]], sep="_")
+      nonlinear_term <- paste(k_term, new_boxes[[1]], sep=" * ")
+      spec[[varname]]$nlt<-nonlinear_term
+      new_diffs[[1]] <- paste(new_diffs[[1]], "-", nonlinear_term)
+      new_parms <- k_term
+      ff <- vector()
+    } 
+
+    # Construct and add FOMC term and add FOMC parameters if needed
+    if(spec[[varname]]$type == "FOMC") {
+      if(match(varname, obs_vars) != 1) {
+        stop("Type FOMC is only possible for the first compartment, which is assumed to be the source compartment")
+      }
+#      if(spec[[varname]]$sink == FALSE) {
+#        stop("Turning off the sink for the FOMC model is not implemented")
+#      }
+      # From p. 53 of the FOCUS kinetics report
+      nonlinear_term <- paste("(alpha/beta) * ((time/beta) + 1)^-1 *", new_boxes[[1]])
+      spec[[varname]]$nlt<-nonlinear_term
+      new_diffs[[1]] <- paste(new_diffs[[1]], "-", nonlinear_term)
+      new_parms <- c("alpha", "beta")
+      ff <- vector()
+    } 
+
+    # Construct and add DFOP term and add DFOP parameters if needed
+    if(spec[[varname]]$type == "DFOP") {
+      if(match(varname, obs_vars) != 1) {
+        stop("Type DFOP is only possible for the first compartment, which is assumed to be the source compartment")
+      }
+#      if(spec[[varname]]$sink == FALSE) {
+#        stop("Turning off the sink for the DFOP model is not implemented")
+#      }
+      # From p. 57 of the FOCUS kinetics report
+      nonlinear_term <- paste("((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 * time)) / (g * exp(-k1 * time) + (1 - g) * exp(-k2 * time))) *", new_boxes[[1]])
+      spec[[varname]]$nlt<-nonlinear_term
+      new_diffs[[1]] <- paste(new_diffs[[1]], "-", nonlinear_term)
+      new_parms <- c("k1", "k2", "g")
+      ff <- vector()
+    } 
+
+    # Construct and add HS term and add HS parameters if needed
+    if(spec[[varname]]$type == "HS") {
+      if(match(varname, obs_vars) != 1) {
+        stop("Type HS is only possible for the first compartment, which is assumed to be the source compartment")
+      }
+#      if(spec[[varname]]$sink == FALSE) {
+#        stop("Turning off the sink for the HS model is not implemented")
+#      }
+      # From p. 55 of the FOCUS kinetics report
+#      nonlinear_term <- paste("ifelse(time <= tb, k1, k2)", "*", new_boxes[[1]])
+      nonlinear_term <- paste("((k1 - k2) / ( 1 + exp( 10*(time - tb) ) ) + k2)", "*", new_boxes[[1]])
+      spec[[varname]]$nlt<-nonlinear_term
+      new_diffs[[1]] <- paste(new_diffs[[1]], "-", nonlinear_term)
+      new_parms <- c("k1", "k2", "tb")
+      ff <- vector()
+    } 
+
+    # Construct terms for transfer to sink and add if appropriate
+
+    #@@@@ REMOVE THIS ?????
+    if(spec[[varname]]$sink) {
+      # Add first-order sink term to first (or only) box for SFO and SFORB
+#      if(spec[[varname]]$type %in% c("SFO", "SFORB")) {
+      if(spec[[varname]]$type == "SFORB") {
+        k_compound_sink <- paste("k", new_boxes[[1]], "sink", sep="_")
+        sink_term <- paste("-", k_compound_sink, "*", new_boxes[[1]])
+        new_diffs[[1]] <- paste(new_diffs[[1]], sink_term)
+        new_parms <- k_compound_sink
+      }
+    }
+   
+    # Add reversible binding if appropriate
+    if(spec[[varname]]$type == "SFORB") {
+      k_free_bound <- paste("k", varname, "free", "bound", sep="_")      
+      k_bound_free <- paste("k", varname, "bound", "free", sep="_")      
+      reversible_binding_terms <- c(
+        paste("-", k_free_bound, "*", new_boxes[[1]], "+", k_bound_free, "*", new_boxes[[2]]),
+        paste("+", k_free_bound, "*", new_boxes[[1]], "-", k_bound_free, "*", new_boxes[[2]]))
+      new_diffs <- paste(new_diffs, reversible_binding_terms)
+      new_parms <- c(new_parms, k_free_bound, k_bound_free)
+    } 
+
+    # Add observed variable to model
+    parms <- c(parms, new_parms)
+    names(new_diffs) <- new_boxes
+    diffs <- c(diffs, new_diffs)
+  }
+  # Transfer between compartments
+  for (varname in obs_vars) {
+    to <- spec[[varname]]$to
+    if(!is.null(to)) {
+      origin_box <- switch(spec[[varname]]$type,
+        SFO = varname,
+        FOMC = varname,
+        DFOP = varname,
+        HS = varname,
+        SFORB = paste(varname, "free", sep="_"))
+      fraction_left <- NULL
+      for (target in to) {
+        target_box <- switch(spec[[target]]$type,
+          SFO = target,
+          SFORB = paste(target, "free", sep="_"))
+        # SFO is no longer special
+        #if(spec[[varname]]$type %in% c("SFO", "SFORB")) {
+        if(spec[[varname]]$type == "SFORB") {
+          k_from_to <- paste("k", origin_box, target_box, sep="_")
+          diffs[[origin_box]] <- paste(diffs[[origin_box]], "-", 
+            k_from_to, "*", origin_box)
+          diffs[[target_box]] <- paste(diffs[[target_box]], "+", 
+            k_from_to, "*", origin_box)
+          parms <- c(parms, k_from_to)
+        }
+        # Handle SFO like the others
+#        if(spec[[varname]]$type %in% c("FOMC", "DFOP", "HS")) {
+        if(spec[[varname]]$type %in% c("FOMC", "DFOP", "HS", "SFO")) {
+          if ( length(to)==1 && !spec[[varname]]$sink ) {
+            # There is only one output, so no need for any flow fractions. Just add the whole flow from the parent
+            diffs[[target_box]] <- paste(diffs[[target_box]], "+", spec[[varname]]$nlt)
+          } else {
+            fraction_to_target = paste("f", varname, "to", target, sep="_")
+            fraction_not_to_target = paste("(1 - ", fraction_to_target, ")", sep="")
+            if(is.null(fraction_left)) {
+              fraction_really_to_target = fraction_to_target
+              fraction_left = fraction_not_to_target
+            } else {
+              # If this is the last output and there is no sink, it gets what's left
+              if ( target==tail(to,1) && !spec[[varname]]$sink ) {
+                 fraction_really_to_target = fraction_left
+              } else {
+# (1-fa)fb version
+#                 fraction_really_to_target = paste(fraction_left, " * ", fraction_to_target, sep="")
+#                 fraction_left = paste(fraction_left, " * ", fraction_not_to_target, sep="")
+# fb version
+                 fraction_really_to_target = fraction_to_target
+                 fraction_left = paste(fraction_left, " - ", fraction_to_target, sep="")
+              }
+            }
+            ff[target_box] = fraction_really_to_target
+            diffs[[target_box]] <- paste(diffs[[target_box]], "+", ff[target_box], "*", spec[[varname]]$nlt)
+            # Add the flow fraction parameter (if it exists)
+            if ( target!=tail(to,1) || spec[[varname]]$sink ) {
+              parms <- c(parms, fraction_to_target)
+            }
+          }
+        }
+      }
+    }
+  }
+  model <- list(diffs = diffs, parms = parms, map = map)
+
+  # Create coefficient matrix if appropriate
+  if (mat) {
+    boxes <- names(diffs)
+    n <- length(boxes)
+    m <- matrix(nrow=n, ncol=n, dimnames=list(boxes, boxes))
+    for (from in boxes) {
+      for (to in boxes) {
+        if (from == to) {
+          #@@@@ OMIT NEXT LINE? !!!!!
+          k.candidate = paste("k", from, c(boxes, "sink"), sep="_")
+          k.candidate = sub("free.*bound", "free_bound", k.candidate)
+          k.candidate = sub("bound.*free", "bound_free", k.candidate)
+          k.effective = intersect(model$parms, k.candidate)
+          m[from,to] = ifelse(length(k.effective) > 0,
+              paste("-", k.effective, collapse = " "), "0")
+        } else {
+          k.candidate = paste("k", from, to, sep="_")
+          k.candidate = sub("free.*bound", "free_bound", k.candidate)
+          k.candidate = sub("bound.*free", "bound_free", k.candidate)
+          k.effective = intersect(model$parms, k.candidate)
+          m[to, from] = ifelse(length(k.effective) > 0,
+              k.effective, "0")
+        }
+      }
+    }
+    model$coefmat <- m
+  }
+
+  if (exists("ff")) model$ff = ff
+  class(model) <- "mkinmod"
+  invisible(model)
+}