From e959fde98f95f3595e01490b67892678bbcd1b27 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Wed, 7 May 2014 14:47:28 +0200 Subject: Fork the gmkin GUI from mkin. See ChangeLog for details --- R/mkinmod.R | 242 ------------------------------------------------------------ 1 file changed, 242 deletions(-) delete mode 100644 R/mkinmod.R (limited to 'R/mkinmod.R') diff --git a/R/mkinmod.R b/R/mkinmod.R deleted file mode 100644 index 2bc6ae3..0000000 --- a/R/mkinmod.R +++ /dev/null @@ -1,242 +0,0 @@ -# Copyright (C) 2010-2013 Johannes Ranke {{{ -# Contact: jranke@uni-bremen.de - -# This file is part of the R package mkin - -# mkin 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 }}} - -mkinmod <- function(..., use_of_ff = "min", speclist = NULL) -{ - if (is.null(speclist)) spec <- list(...) - else spec <- speclist - obs_vars <- names(spec) - - # Check if any of the names of the observed variables contains any other - for (obs_var in obs_vars) { - if (length(grep(obs_var, obs_vars)) > 1) stop("Sorry, variable names can not contain each other") - if (grepl("_to_", obs_var)) stop("Sorry, names of observed variables can not contain _to_") - } - - if (!use_of_ff %in% c("min", "max")) - stop("The use of formation fractions 'use_of_ff' can only be 'min' or 'max'") - - # 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() - # }}} - - # Give a warning when a model with time dependent degradation uses formation {{{ - # fractions - if(spec[[1]]$type %in% c("FOMC", "DFOP", "HS")) { - mat = FALSE - } else mat = TRUE - #}}} - - # Establish a list of differential equations as well as a map from observed {{{ - # compartments to differential equations - for (varname in obs_vars) - { - # Check the type component of the compartment specification {{{ - if(is.null(spec[[varname]]$type)) stop( - "Every part of the model specification 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") - 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, - 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="") - names(new_diffs) <- new_boxes - diffs <- c(diffs, new_diffs) #}}} - } #}}} - - # Create content of differential equations and build parameter list {{{ - for (varname in obs_vars) - { - # 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 this is not explicitly excluded by the user by - # specifying sink=FALSE - if(is.null(spec[[varname]]$sink)) spec[[varname]]$sink <- TRUE - if(spec[[varname]]$type %in% c("SFO", "SFORB")) { # {{{ Add SFO or SFORB decline term - if (use_of_ff == "min") { # Minimum use of formation fractions - if(spec[[varname]]$sink) { - # If sink is required, add first-order sink term - k_compound_sink <- paste("k", box_1, "sink", sep="_") - parms <- c(parms, k_compound_sink) - decline_term <- paste(k_compound_sink, "*", box_1) - } else { # otherwise no decline term needed here - decline_term = "0" - } - } else { - k_compound <- paste("k", box_1, sep="_") - parms <- c(parms, k_compound) - decline_term <- paste(k_compound, "*", box_1) - } - } #}}} - if(spec[[varname]]$type == "FOMC") { # {{{ Add FOMC decline term - # From p. 53 of the FOCUS kinetics report - decline_term <- paste("(alpha/beta) * ((time/beta) + 1)^-1 *", box_1) - parms <- c(parms, "alpha", "beta") - } #}}} - if(spec[[varname]]$type == "DFOP") { # {{{ Add DFOP decline term - # From p. 57 of the FOCUS kinetics report - decline_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") - } #}}} - if(spec[[varname]]$type == "HS") { # {{{ Add HS decline term - # From p. 55 of the FOCUS kinetics report - decline_term <- paste("ifelse(time <= tb, k1, k2)", "*", box_1) - parms <- c(parms, "k1", "k2", "tb") - } #}}} - # Add origin decline term to box 1 (usually the only box, unless type is SFORB)#{{{ - diffs[[box_1]] <- paste(diffs[[box_1]], "-", decline_term)#}}} - if(spec[[varname]]$type == "SFORB") { # {{{ Add SFORB reversible binding terms - box_2 = map[[varname]][[2]] - if (use_of_ff == "min") { # Minimum use of formation fractions - k_free_bound <- paste("k", varname, "free", "bound", sep="_") - k_bound_free <- paste("k", varname, "bound", "free", sep="_") - parms <- c(parms, k_free_bound, k_bound_free) - reversible_binding_term_1 <- paste("-", k_free_bound, "*", box_1, "+", - k_bound_free, "*", box_2) - reversible_binding_term_2 <- paste("+", k_free_bound, "*", box_1, "-", - k_bound_free, "*", box_2) - } else { # Use formation fractions also for the free compartment - stop("The maximum use of formation fractions is not supported for SFORB models") - # The problems were: Calculation of dissipation times did not work in this case - # and the coefficient matrix is not generated correctly by the code present - # in this file in this case - f_free_bound <- paste("f", varname, "free", "bound", sep="_") - k_bound_free <- paste("k", varname, "bound", "free", sep="_") - parms <- c(parms, f_free_bound, k_bound_free) - reversible_binding_term_1 <- paste("+", k_bound_free, "*", box_2) - reversible_binding_term_2 <- paste("+", f_free_bound, "*", k_compound, "*", box_1, "-", - k_bound_free, "*", box_2) - } - diffs[[box_1]] <- paste(diffs[[box_1]], reversible_binding_term_1) - diffs[[box_2]] <- paste(diffs[[box_2]], reversible_binding_term_2) - } #}}} - - # Transfer between compartments#{{{ - to <- spec[[varname]]$to - if(!is.null(to)) { - # 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 (use_of_ff == "min" && spec[[varname]]$type %in% c("SFO", - "SFORB")) { - k_from_to <- paste("k", origin_box, target_box, sep="_") - parms <- c(parms, k_from_to) - diffs[[origin_box]] <- paste(diffs[[origin_box]], "-", - k_from_to, "*", origin_box) - diffs[[target_box]] <- paste(diffs[[target_box]], "+", - k_from_to, "*", origin_box) - } else { - if (!spec[[varname]]$sink) { - stop("Turning off the sink when using formation fractions is not supported") - } - fraction_to_target = paste("f", origin_box, "to", target, sep="_") - parms <- c(parms, fraction_to_target) - diffs[[target_box]] <- paste(diffs[[target_box]], "+", - fraction_to_target, "*", decline_term) - } - } - } #}}} - } #}}} - - model <- list(diffs = diffs, parms = parms, map = map, spec = spec, use_of_ff = use_of_ff) - - # Create coefficient matrix if appropriate#{{{ - if (mat) { - boxes <- names(diffs) - n <- length(boxes) - m <- matrix(nrow=n, ncol=n, dimnames=list(boxes, boxes)) - - 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 - }#}}} - - class(model) <- "mkinmod" - return(model) -} -# vim: set foldmethod=marker ts=2 sw=2 expandtab: -- cgit v1.2.1