Start of the transition to fitting transformed parameters.

Many things are broken now (see TODO list)



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

1 files changed, 61 insertions, 107 deletions

diff --git a/R/mkinmod.R b/R/mkinmod.R
index 31c778cb..b49c813b 100644
--- a/R/mkinmod.R
+++ b/R/mkinmod.R
@@ -1,7 +1,7 @@
 # $Id: mkinmod.R 71 2010-09-12 01:13:36Z jranke $

 

-# Copyright (C) 2010 Johannes Ranke

-# Contact: mkin-devel@lists.berlios.de

+# Copyright (C) 2010-2012 Johannes Ranke

+# Contact: jranke@uni-bremen.de

 

 # This file is part of the R package mkin

 

@@ -33,21 +33,27 @@ mkinmod <- function(...)
   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",

+      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

 

   # Establish list of differential equations

+  # as well as map from observed compartments to 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()

+    if(spec[[varname]]$type %in% c("FOMC", "DFOP", "HS") & match(varname, obs_vars) != 1) {

+        stop(paste("Types FOMC, DFOP and HS are only implemented for the first compartment,", 

+                   "which is assumed to be the source compartment"))

+    }

 

     # New (sub)compartments (boxes) needed for the model type

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

@@ -62,124 +68,75 @@ mkinmod <- function(...)
 

     # Start a new differential equation for each new box

     new_diffs <- paste("d_", new_boxes, " =", sep="")

+    names(new_diffs) <- new_boxes

+    diffs <- c(diffs, new_diffs)

+  }

+

+  # Create content of differential equations and build parameter list

+  for (varname in obs_vars)

+  {

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

+    box_1 = map[[varname]][[1]] # This is the only box unless type is SFORB

+    if(spec[[varname]]$type %in% c("SFO", "SFORB")) {

+      k_compound <- paste("k", box_1, sep="_")

+      origin_term <- paste(k_compound, "*", box_1)

+      parms <- c(parms, k_compound)

+    }

 

-    # Turn on sink if not specified otherwise

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

 

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

-      new_diffs[[1]] <- paste(new_diffs[[1]], "-", nonlinear_term)

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

-      ff <- vector()

+      origin_term <- paste("(alpha/beta) * ((time/beta) + 1)^-1 *", box_1)

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

     } 

 

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

-      new_diffs[[1]] <- paste(new_diffs[[1]], "-", nonlinear_term)

-      new_parms <- c("k1", "k2", "g")

-      ff <- vector()

+      origin_term <- paste("((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 * time)) / (g * exp(-k1 * time) + (1 - g) * exp(-k2 * time))) *", box_1)

+      parms <- c(parms, "k1", "k2", "g")

     } 

 

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

-      new_diffs[[1]] <- paste(new_diffs[[1]], "-", nonlinear_term)

-      new_parms <- c("k1", "k2", "tb")

-      ff <- vector()

+      origin_term <- paste("ifelse(time <= tb, k1, k2)", "*", box_1)

+      parms <- c(parms, "k1", "k2", "tb")

     } 

 

-    # Construct terms for transfer to sink and add if appropriate

-

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

-        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 origin decline term to box 1 (usually the only box, unless type is SFORB)

+    diffs[[box_1]] <- paste(diffs[[box_1]], "-", origin_term)

    

-    # Add reversible binding if appropriate

+    # Add reversible binding for SFORB models

     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)

+      box_2 = map[[varname]][[2]]

+      k_free_bound <- paste("k", varname, "free", "bound", sep="_")

+      k_bound_free <- paste("k", varname, "bound", "free", sep="_")

+      reversible_binding_term_1 <- paste("-", k_free_bound, "*", box_1, "+",

+        k_bound_free, "*", box_2)

+      diffs[[box_1]] <- paste(diffs[[box_1]], reversible_binding_term_1)

+      reversible_binding_term_2 <- paste("+", k_free_bound, "*", box_1, "-",

+        k_bound_free, "*", box_2)

+      diffs[[box_2]] <- paste(diffs[[box_2]], reversible_binding_term_2)

+      parms <- c(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) {

+    # Transfer between compartments

     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

+      # Name of box from which transfer takes place

+      origin_box <- box_1

+

+      # Add transfer terms to listed compartments

       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 %in% c("FOMC", "DFOP", "HS")) {

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

-          }

-          ff[target_box] = fraction_really_to_target

-          diffs[[target_box]] <- paste(diffs[[target_box]], "+", 

-            ff[target_box], "*", nonlinear_term)

-          parms <- c(parms, fraction_to_target)

-        }

+        fraction_to_target = paste("f", origin_box, "to", target, sep="_")

+        diffs[[target_box]] <- paste(diffs[[target_box]], "+", 

+            fraction_to_target, "*", origin_term)

+        parms <- c(parms, fraction_to_target)

       }

     }

   }

@@ -193,26 +150,23 @@ mkinmod <- function(...)
     for (from in boxes) {

       for (to in boxes) {

         if (from == to) {

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

+          k.candidate = paste("k", from, sep="_")

+          m[from,to] = ifelse(k.candidate %in% model$parms,

+              paste("-", k.candidate), "0")

         } else {

+          f.candidate = paste("f", from, "to", to, sep="_")

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

+          m[to, from] = ifelse(f.candidate %in% model$parms,

+            paste(f.candidate, " * k_", from, sep=""), 

+            ifelse(k.candidate %in% model$parms, k.candidate, "0"))

         }

       }

     }

     model$coefmat <- m

   }

 

-  if (exists("ff")) model$ff = ff

   class(model) <- "mkinmod"

-  invisible(model)

+  return(model)

 }