Version 2.0v2.0

1 files changed, 75 insertions, 25 deletions

diff --git a/CakeModel.R b/CakeModel.R
index a93d7af..ce6bc01 100644
--- a/CakeModel.R
+++ b/CakeModel.R
@@ -1,5 +1,5 @@
 # Was Id: mkinmod.R 71 2010-09-12 01:13:36Z jranke 
-# $Id: CakeModel.R 216 2011-07-05 14:35:03Z nelr $
+# $Id$
 
 # The CAKE R modules are based on mkin
 # Portions Johannes Ranke, 2010
@@ -22,12 +22,13 @@
 #    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(...)
+CakeModel <- function(..., use_of_ff = "max")
 {
   spec <- list(...)
   obs_vars <- names(spec)
 
-  # The returned model will be a list of character vectors, containing
+  if (!use_of_ff %in% c("min", "max"))
+    stop("The use of formation fractions 'use_of_ff' can only be 'min' or 'max'")
   # differential equations, parameter names and a mapping from model variables
   # to observed variables. If possible, a matrix representation of the 
   # differential equations is included
@@ -69,21 +70,39 @@ CakeModel <- function(...)
 
     # Start a new differential equation for each new box
     new_diffs <- paste("d_", new_boxes, " =", sep="")
-
+    
+    # Get the name of the box(es) we are working on for the decline term(s)
+    box_1 = map[[varname]][[1]] # This is the only box unless type is SFORB
+    
     # 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") {
+      if (use_of_ff == "min") { # Minimum use of formation fractions
+        if(spec[[varname]]$sink) {
       # From p. 53 of the FOCUS kinetics report
+      k_term <- paste("k", new_boxes[[1]], "sink", 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()
+      decline_term <- paste(nonlinear_term, "*", new_boxes[[1]])
+    } else { # otherwise no decline term needed here
+      decline_term = "0" 
+    }
+  } else {
       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()
-    } 
+    decline_term <- paste(nonlinear_term, "*", new_boxes[[1]])
+  }
+ } 
 
     # Construct and add FOMC term and add FOMC parameters if needed
     if(spec[[varname]]$type == "FOMC") {
@@ -225,33 +244,64 @@ CakeModel <- function(...)
       }
     }
   }
-  model <- list(diffs = diffs, parms = parms, map = map)
+  model <- list(diffs = diffs, parms = parms, map = map, use_of_ff = use_of_ff)
 
-  # Create coefficient matrix if appropriate
+  # 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")
+
+    if (use_of_ff == "min") { # Minimum use of formation fractions
+      for (from in boxes) {
+        for (to in boxes) {
+          if (from == to) { # diagonal elements
+            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 {          # off-diagonal elements
+            k.candidate = paste("k", from, to, sep="_")
+	    if (sub("_free$", "", from) == sub("_bound$", "", to)) {
+              k.candidate = paste("k", sub("_free$", "_free_bound", from), sep="_")
+	    }
+	    if (sub("_bound$", "", from) == sub("_free$", "", to)) {
+              k.candidate = paste("k", sub("_bound$", "_bound_free", from), sep="_")
+	    }
+            k.effective = intersect(model$parms, k.candidate)
+            m[to, from] = ifelse(length(k.effective) > 0,
+                k.effective, "0")
+          }
+        }
+      } 
+    } else {  # Use formation fractions where possible
+      for (from in boxes) {
+        for (to in boxes) {
+          if (from == to) { # diagonal elements
+            k.candidate = paste("k", from, sep="_")
+            m[from,to] = ifelse(k.candidate %in% model$parms,
+                paste("-", k.candidate), "0")
+            if(grepl("_free", from)) { # add transfer to bound compartment for SFORB
+              m[from,to] = paste(m[from,to], "-", paste("k", from, "bound", sep="_"))
+            }
+            if(grepl("_bound", from)) { # add backtransfer to free compartment for SFORB
+              m[from,to] = paste("- k", from, "free", sep="_")
+            }
+            m[from,to] = m[from,to]
+          } else {          # off-diagonal elements
+            f.candidate = paste("f", from, "to", to, sep="_")
+            k.candidate = paste("k", from, to, sep="_")
+	    # SFORB with maximum use of formation fractions not implemented, see above
+            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
   }