Much more complete version that was just submitted to CRAN.

git-svn-id: svn+ssh://svn.r-forge.r-project.org/svnroot/kinfit/pkg/mkin@9 edb9625f-4e0d-4859-8d74-9fd3b1da38cb

1 files changed, 49 insertions, 17 deletions

diff --git a/R/mkinmod.R b/R/mkinmod.R
index 7e64f879..643e855e 100644
--- a/R/mkinmod.R
+++ b/R/mkinmod.R
@@ -1,38 +1,53 @@
-mkinmod <- function(spec = list(parent = list(type = "SFO", to = NA, sink = TRUE)))

+mkinmod <- function(...)

 {

-  if(!is.list(spec)) stop("spec must be a list of model specifications for each observed variable")

+  spec <- list(...)

+  obs_vars <- names(spec)

 

-  # The returned model will be a list of two character vectors, containing

-  # differential equations and parameter names

+  # The returned model will be a list of character vectors, containing

+  # differential equations, parameter names and a mapping from model variables

+  # to observed variables

   parms <- vector()

   diffs <- vector()

   map <- list()

 

   # Establish list of differential equations

-  for (varname in names(spec))

+  for (varname in obs_vars)

   {

-    if(!spec[[varname]]$type %in% c("SFO", "SFORB")) stop("Available types are SFO and SFORB only")

+    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", "SFORB")) stop("Available types are SFO, FOMC 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,

       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="")

 

     # Construct terms for transfer to sink and add if appropriate

+    if(is.null(spec[[varname]]$sink)) spec[[varname]]$sink <- TRUE

     if(spec[[varname]]$sink) {

-      k_compound_sink <- paste("k", new_boxes[[1]], "sink", sep="_")

-      sink_term <- paste("-", k_compound_sink, "*", new_boxes[[1]])

-      # Only add sink term to first (or only) box for SFO and SFORB models

+      # Add first-order sink term to first (or only) box for SFO and SFORB models

       if(spec[[varname]]$type %in% c("SFO", "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

+      }

+      if(spec[[varname]]$type == "FOMC") {

+        if(match(varname, obs_vars) != 1) {

+          stop("Type FOMC is only allowed for the first compartment, which is assumed to be the source compartment")

+        }

+        # From p. 53 of the FOCUS kinetics report

+        fomc_term <- paste("(alpha/beta) * ((time/beta) + 1)^-1 *", new_boxes[[1]])

+        new_diffs[[1]] <- paste(new_diffs[[1]], "-", fomc_term)

+        new_parms <- c("alpha", "beta")

       }

-      new_parms <- k_compound_sink

     }

    

     # Add reversible binding if appropriate

@@ -53,24 +68,41 @@ mkinmod <- function(spec = list(parent = list(type = "SFO", to = NA, sink = TRUE
   }

 

   # Transfer between compartments

-  for (varname in names(spec)) {

+  for (varname in obs_vars) {

     to <- spec[[varname]]$to

-    if(!is.na(to)) {

+    if(!is.null(to)) {

       origin_box <- switch(spec[[varname]]$type,

         SFO = varname,

+        FOMC = 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="_"))

+        if(spec[[varname]]$type %in% c("SFO", "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)

+        }

+        if(spec[[varname]]$type == "FOMC") {

+          fraction_to_target = paste("f_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 {

+            fraction_really_to_target = paste(fraction_left, " * ", fraction_to_target, sep="")

+            fraction_left = paste(fraction_left, " * ", fraction_not_to_target, sep="")

+          }

+          diffs[[target_box]] <- paste(diffs[[target_box]], "+", fraction_really_to_target, "*", fomc_term)

+          parms <- c(parms, fraction_to_target)

+        }

       }

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

     }

   }

   model <- list(diffs = diffs, parms = parms, map = map)

   class(model) <- "mkinmod"

-  return(model)

+  invisible(model)

 }