diff options
Diffstat (limited to 'CakeModel.R')
| -rw-r--r-- | CakeModel.R | 100 |
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 } |