Use roxygen for functions and methods

1 files changed, 483 insertions, 383 deletions

diff --git a/R/mkinmod.R b/R/mkinmod.R
index 26148f18..a1ae0021 100644
--- a/R/mkinmod.R
+++ b/R/mkinmod.R
@@ -1,383 +1,483 @@
-# Copyright (C) 2010-2015,2019 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 <http://www.gnu.org/licenses/> }}}

-

-mkinmod <- function(..., use_of_ff = "min", speclist = NULL, quiet = FALSE, verbose = FALSE)

-{

-  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 (obs_var == "sink") stop("Naming a compound 'sink' is not supported")

-  }

-

-  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 (if supported), parameter names and a mapping from

-  # model variables to observed variables. If possible, a matrix representation

-  # of the differential equations is included

-  # Compiling the functions from the C code generated below only works if the

-  # implicit assumption about differential equations specified below

-  # is satisfied

-  parms <- vector()

-  # }}}

-

-  # Do not return a coefficient matrix mat when FOMC, IORE, DFOP, HS or logistic is used for the parent {{{

-  if(spec[[1]]$type %in% c("FOMC", "IORE", "DFOP", "HS", "logistic")) {

-    mat = FALSE

-  } else mat = TRUE

-  #}}}

-

-  # Establish a list of differential equations as well as a map from observed {{{

-  # compartments to differential equations

-  diffs <- vector()

-  map <- list()

-  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", "IORE", "DFOP", "HS", "SFORB", "logistic")) stop(

-      "Available types are SFO, FOMC, IORE, DFOP, HS, SFORB and logistic only")

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

-        stop(paste("Types FOMC, DFOP, HS and logistic 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,

-      IORE = varname,

-      DFOP = varname,

-      HS = varname,

-      logistic = 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", "IORE", "SFORB")) { # {{{ Add decline term

-      if (use_of_ff == "min") { # Minimum use of formation fractions

-        if(spec[[varname]]$type == "IORE" && length(spec[[varname]]$to) > 0) {

-           stop("Transformation reactions from compounds modelled with IORE\n",

-                "are only supported with formation fractions (use_of_ff = 'max')")

-        }

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

-          # If sink is required, add first-order/IORE sink term

-          k_compound_sink <- paste("k", box_1, "sink", sep = "_")

-          if(spec[[varname]]$type == "IORE") {

-            k_compound_sink <- paste("k__iore", box_1, "sink", sep = "_")

-          }

-          parms <- c(parms, k_compound_sink)

-          decline_term <- paste(k_compound_sink, "*", box_1)

-          if(spec[[varname]]$type == "IORE") {

-            N <- paste("N", box_1, sep = "_")

-            parms <- c(parms, N)

-            decline_term <- paste0(decline_term, "^", N)

-          }

-        } else { # otherwise no decline term needed here

-          decline_term = "0"

-        }

-      } else {

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

-        if(spec[[varname]]$type == "IORE") {

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

-        }

-        parms <- c(parms, k_compound)

-        decline_term <- paste(k_compound, "*", box_1)

-        if(spec[[varname]]$type == "IORE") {

-          N <- paste("N", box_1, sep = "_")

-          parms <- c(parms, N)

-          decline_term <- paste0(decline_term, "^", N)

-        }

-      }

-    } #}}}

-    if(spec[[varname]]$type == "FOMC") { # {{{ Add FOMC decline term

-      # From p. 53 of the FOCUS kinetics report, without the power function so it works in C

-      decline_term <- paste("(alpha/beta) * 1/((time/beta) + 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")

-    } #}}}

-    HS_decline <- "ifelse(time <= tb, k1, k2)" # Used below for automatic translation to C

-    if(spec[[varname]]$type == "HS") { # {{{ Add HS decline term

-      # From p. 55 of the FOCUS kinetics report

-      decline_term <- paste(HS_decline, "*", box_1)

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

-    } #}}}

-    if(spec[[varname]]$type == "logistic") { # {{{ Add logistic decline term

-      # From p. 67 of the FOCUS kinetics report (2014)

-      decline_term <- paste("(k0 * kmax)/(k0 + (kmax - k0) * exp(-r * time)) *", box_1)

-      parms <- c(parms, "kmax", "k0", "r")

-    } #}}}

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

-

-      # Number of targets

-      n_targets = length(to)

-

-      # Add transfer terms to listed compartments

-      for (target in to) {

-        if (!target %in% obs_vars) stop("You did not specify a submodel for target variable ", target)

-        target_box <- switch(spec[[target]]$type,

-          SFO = target,

-          IORE = 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 {

-          # Do not introduce a formation fraction if this is the only target

-          if (spec[[origin_box]]$sink == FALSE && n_targets == 1) {

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

-                                         decline_term)

-          } else {

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

-            # Special case: singular pathway and no sink

-            if (spec[[from]]$sink == FALSE && length(spec[[from]]$to) == 1 && to %in% spec[[from]]$to) {

-              m[to, from] = paste("k", from, sep = "_")

-            }

-          }

-        }

-      }

-    } # }}}

-    model$coefmat <- m

-  }#}}}

-

-  # Try to create a function compiled from C code if more than one observed {{{

-  # variable and gcc is available

-  if (length(obs_vars) > 1) {

-    if (Sys.which("gcc") != "") {

-

-      # Translate the R code for the derivatives to C code

-      diffs.C <- paste(diffs, collapse = ";\n")

-      diffs.C <- paste0(diffs.C, ";")

-

-      # HS

-      diffs.C <- gsub(HS_decline, "(time <= tb ? k1 : k2)", diffs.C, fixed = TRUE)

-

-      for (i in seq_along(diffs)) {

-        state_var <- names(diffs)[i]

-

-        # IORE

-        if (state_var %in% obs_vars) {

-          if (spec[[state_var]]$type == "IORE") {

-            diffs.C <- gsub(paste0(state_var, "^N_", state_var),

-                            paste0("pow(y[", i - 1, "], N_", state_var, ")"),

-                            diffs.C, fixed = TRUE)

-          }

-        }

-

-        # Replace d_... terms by f[i-1]

-        # First line

-        pattern <- paste0("^d_", state_var)

-        replacement <- paste0("\nf[", i - 1, "]")

-        diffs.C <- gsub(pattern, replacement, diffs.C)

-        # Other lines

-        pattern <- paste0("\\nd_", state_var)

-        replacement <- paste0("\nf[", i - 1, "]")

-        diffs.C <- gsub(pattern, replacement, diffs.C)

-

-        # Replace names of observed variables by y[i],

-        # making the implicit assumption that the observed variables only occur after "* "

-        pattern <- paste0("\\* ", state_var)

-        replacement <- paste0("* y[", i - 1, "]")

-        diffs.C <- gsub(pattern, replacement, diffs.C)

-      }

-

-      derivs_sig <- signature(n = "integer", t = "numeric", y = "numeric",

-                              f = "numeric", rpar = "numeric", ipar = "integer")

-

-      # Declare the time variable in the body of the function if it is used

-      derivs_code <- if (spec[[1]]$type %in% c("FOMC", "DFOP", "HS")) {

-        paste0("double time = *t;\n", diffs.C)

-      } else {

-        diffs.C

-      }

-

-      # Define the function initializing the parameters

-      npar <- length(parms)

-      initpar_code <- paste0(

-        "static double parms [", npar, "];\n",

-        paste0("#define ", parms, " parms[", 0:(npar - 1), "]\n", collapse = ""),

-        "\n",

-        "void initpar(void (* odeparms)(int *, double *)) {\n",

-        "    int N = ", npar, ";\n",

-        "    odeparms(&N, parms);\n",

-        "}\n\n")

-

-      # Try to build a shared library

-      cf <- try(cfunction(list(func = derivs_sig), derivs_code,

-                          otherdefs = initpar_code,

-                          verbose = verbose,

-                          convention = ".C", language = "C"),

-                silent = TRUE)

-

-      if (!inherits(cf, "try-error")) {

-        if (!quiet) message("Successfully compiled differential equation model from auto-generated C code.")

-        model$cf <- cf

-      }

-    }

-  }

-  # }}}

-

-  class(model) <- "mkinmod"

-  return(model)

-}

-

-print.mkinmod <- function(x, ...) {

-  cat("<mkinmod> model generated with\n")

-  cat("Use of formation fractions $use_of_ff:", x$use_of_ff, "\n")

-  cat("Specification $spec:\n")

-  for (obs in names(x$spec)) {

-    cat("$", obs, "\n", sep = "")

-    spl <- x$spec[[obs]]

-    cat("$type:", spl$type)

-    if (!is.null(spl$to) && length(spl$to)) cat("; $to: ", paste(spl$to, collapse = ", "), sep = "")

-    cat("; $sink: ", spl$sink, sep = "")

-    if (!is.null(spl$full_name)) if (!is.na(spl$full_name)) cat("; $full_name:", spl$full_name)

-    cat("\n")

-  }

-  if (is.matrix(x$coefmat)) cat("Coefficient matrix $coefmat available\n")

-  if (!is.null(x$cf)) cat("Compiled model $cf available\n")

-  cat("Differential equations:\n")

-  nice_diffs <- gsub("^(d.*) =", "\\1/dt =", x[["diffs"]])

-  writeLines(strwrap(nice_diffs, exdent = 11))

-}

-# vim: set foldmethod=marker ts=2 sw=2 expandtab:

+#' Function to set up a kinetic model with one or more state variables
+#'
+#' The function usually takes several expressions, each assigning a compound
+#' name to a list, specifying the kinetic model type and reaction or transfer
+#' to other observed compartments. Instead of specifying several expressions, a
+#' list of lists can be given in the speclist argument.
+#'
+#' For the definition of model types and their parameters, the equations given
+#' in the FOCUS and NAFTA guidance documents are used.
+#'
+#' @param ...  For each observed variable, a list has to be specified as an
+#'   argument, containing at least a component \code{type}, specifying the type
+#'   of kinetics to use for the variable. Currently, single first order
+#'   kinetics "SFO", indeterminate order rate equation kinetics "IORE", or
+#'   single first order with reversible binding "SFORB" are implemented for all
+#'   variables, while "FOMC", "DFOP" and "HS" can additionally be chosen for
+#'   the first variable which is assumed to be the source compartment.
+#'   Additionally, each component of the list can include a character vector
+#'   \code{to}, specifying names of variables to which a transfer is to be
+#'   assumed in the model.  If the argument \code{use_of_ff} is set to "min"
+#'   (default) and the model for the compartment is "SFO" or "SFORB", an
+#'   additional component of the list can be "sink=FALSE" effectively fixing
+#'   the flux to sink to zero.
+#' @param use_of_ff Specification of the use of formation fractions in the
+#'   model equations and, if applicable, the coefficient matrix. If "min", a
+#'   minimum use of formation fractions is made in order to avoid fitting the
+#'   product of formation fractions and rate constants. If "max", formation
+#'   fractions are always used.
+#' @param speclist The specification of the observed variables and their
+#'   submodel types and pathways can be given as a single list using this
+#'   argument. Default is NULL.
+#' @param quiet Should messages be suppressed?
+#' @param verbose If \code{TRUE}, passed to \code{\link{cfunction}} if
+#'   applicable to give detailed information about the C function being built.
+#' @importFrom methods signature
+#' @importFrom inline cfunction
+#' @return A list of class \code{mkinmod} for use with \code{\link{mkinfit}},
+#'   containing, among others,
+#'   \item{diffs}{
+#'     A vector of string representations of differential equations, one for
+#'     each modelling variable.
+#'   }
+#'   \item{map}{
+#'     A list containing named character vectors for each observed variable,
+#'     specifying the modelling variables by which it is represented.
+#'   }
+#'   \item{use_of_ff}{
+#'     The content of \code{use_of_ff} is passed on in this list component.
+#'   }
+#'   \item{coefmat}{
+#'     The coefficient matrix, if the system of differential equations can be
+#'     represented by one.
+#'   }
+#'   \item{ll}{
+#'     The likelihood function, taking the parameter vector as the first argument.
+#'   }
+#' @note The IORE submodel is not well tested for metabolites. When using this
+#'   model for metabolites, you may want to read the second note in the help
+#'   page to \code{\link{mkinfit}}.
+#' @author Johannes Ranke
+#' @references FOCUS (2006) \dQuote{Guidance Document on Estimating Persistence
+#'   and Degradation Kinetics from Environmental Fate Studies on Pesticides in
+#'   EU Registration} Report of the FOCUS Work Group on Degradation Kinetics,
+#'   EC Document Reference Sanco/10058/2005 version 2.0, 434 pp,
+#'   \url{http://esdac.jrc.ec.europa.eu/projects/degradation-kinetics}
+#'
+#'   NAFTA Technical Working Group on Pesticides (not dated) Guidance for
+#'   Evaluating and Calculating Degradation Kinetics in Environmental Media
+#' @examples
+#'
+#' # Specify the SFO model (this is not needed any more, as we can now mkinfit("SFO", ...)
+#' SFO <- mkinmod(parent = list(type = "SFO"))
+#'
+#' # One parent compound, one metabolite, both single first order
+#' SFO_SFO <- mkinmod(
+#'   parent = mkinsub("SFO", "m1"),
+#'   m1 = mkinsub("SFO"))
+#'
+#' \dontrun{
+#' # The above model used to be specified like this, before the advent of mkinsub()
+#' SFO_SFO <- mkinmod(
+#'   parent = list(type = "SFO", to = "m1"),
+#'   m1 = list(type = "SFO"))
+#'
+#' # Show details of creating the C function
+#' SFO_SFO <- mkinmod(
+#'   parent = mkinsub("SFO", "m1"),
+#'   m1 = mkinsub("SFO"), verbose = TRUE)
+#'
+#' # If we have several parallel metabolites
+#' # (compare tests/testthat/test_synthetic_data_for_UBA_2014.R)
+#' m_synth_DFOP_par <- mkinmod(parent = mkinsub("DFOP", c("M1", "M2")),
+#'                            M1 = mkinsub("SFO"),
+#'                            M2 = mkinsub("SFO"),
+#'                            use_of_ff = "max", quiet = TRUE)
+#'
+#' fit_DFOP_par_c <- mkinfit(m_synth_DFOP_par,
+#'                           synthetic_data_for_UBA_2014[[12]]$data,
+#'                           quiet = TRUE)
+#' }
+#'
+#' @export mkinmod
+mkinmod <- function(..., use_of_ff = "min", speclist = NULL, quiet = FALSE, verbose = FALSE)
+{
+  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 (obs_var == "sink") stop("Naming a compound 'sink' is not supported")
+  }
+
+  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 (if supported), parameter names and a mapping from
+  # model variables to observed variables. If possible, a matrix representation
+  # of the differential equations is included
+  # Compiling the functions from the C code generated below only works if the
+  # implicit assumption about differential equations specified below
+  # is satisfied
+  parms <- vector()
+  # }}}
+
+  # Do not return a coefficient matrix mat when FOMC, IORE, DFOP, HS or logistic is used for the parent {{{
+  if(spec[[1]]$type %in% c("FOMC", "IORE", "DFOP", "HS", "logistic")) {
+    mat = FALSE
+  } else mat = TRUE
+  #}}}
+
+  # Establish a list of differential equations as well as a map from observed {{{
+  # compartments to differential equations
+  diffs <- vector()
+  map <- list()
+  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", "IORE", "DFOP", "HS", "SFORB", "logistic")) stop(
+      "Available types are SFO, FOMC, IORE, DFOP, HS, SFORB and logistic only")
+    if(spec[[varname]]$type %in% c("FOMC", "DFOP", "HS", "logistic") & match(varname, obs_vars) != 1) {
+        stop(paste("Types FOMC, DFOP, HS and logistic 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,
+      IORE = varname,
+      DFOP = varname,
+      HS = varname,
+      logistic = 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", "IORE", "SFORB")) { # {{{ Add decline term
+      if (use_of_ff == "min") { # Minimum use of formation fractions
+        if(spec[[varname]]$type == "IORE" && length(spec[[varname]]$to) > 0) {
+           stop("Transformation reactions from compounds modelled with IORE\n",
+                "are only supported with formation fractions (use_of_ff = 'max')")
+        }
+        if(spec[[varname]]$sink) {
+          # If sink is required, add first-order/IORE sink term
+          k_compound_sink <- paste("k", box_1, "sink", sep = "_")
+          if(spec[[varname]]$type == "IORE") {
+            k_compound_sink <- paste("k__iore", box_1, "sink", sep = "_")
+          }
+          parms <- c(parms, k_compound_sink)
+          decline_term <- paste(k_compound_sink, "*", box_1)
+          if(spec[[varname]]$type == "IORE") {
+            N <- paste("N", box_1, sep = "_")
+            parms <- c(parms, N)
+            decline_term <- paste0(decline_term, "^", N)
+          }
+        } else { # otherwise no decline term needed here
+          decline_term = "0"
+        }
+      } else {
+        k_compound <- paste("k", box_1, sep = "_")
+        if(spec[[varname]]$type == "IORE") {
+          k_compound <- paste("k__iore", box_1, sep = "_")
+        }
+        parms <- c(parms, k_compound)
+        decline_term <- paste(k_compound, "*", box_1)
+        if(spec[[varname]]$type == "IORE") {
+          N <- paste("N", box_1, sep = "_")
+          parms <- c(parms, N)
+          decline_term <- paste0(decline_term, "^", N)
+        }
+      }
+    } #}}}
+    if(spec[[varname]]$type == "FOMC") { # {{{ Add FOMC decline term
+      # From p. 53 of the FOCUS kinetics report, without the power function so it works in C
+      decline_term <- paste("(alpha/beta) * 1/((time/beta) + 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")
+    } #}}}
+    HS_decline <- "ifelse(time <= tb, k1, k2)" # Used below for automatic translation to C
+    if(spec[[varname]]$type == "HS") { # {{{ Add HS decline term
+      # From p. 55 of the FOCUS kinetics report
+      decline_term <- paste(HS_decline, "*", box_1)
+      parms <- c(parms, "k1", "k2", "tb")
+    } #}}}
+    if(spec[[varname]]$type == "logistic") { # {{{ Add logistic decline term
+      # From p. 67 of the FOCUS kinetics report (2014)
+      decline_term <- paste("(k0 * kmax)/(k0 + (kmax - k0) * exp(-r * time)) *", box_1)
+      parms <- c(parms, "kmax", "k0", "r")
+    } #}}}
+    # 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
+
+      # Number of targets
+      n_targets = length(to)
+
+      # Add transfer terms to listed compartments
+      for (target in to) {
+        if (!target %in% obs_vars) stop("You did not specify a submodel for target variable ", target)
+        target_box <- switch(spec[[target]]$type,
+          SFO = target,
+          IORE = 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 {
+          # Do not introduce a formation fraction if this is the only target
+          if (spec[[origin_box]]$sink == FALSE && n_targets == 1) {
+            diffs[[target_box]] <- paste(diffs[[target_box]], "+",
+                                         decline_term)
+          } else {
+            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"))
+            # Special case: singular pathway and no sink
+            if (spec[[from]]$sink == FALSE && length(spec[[from]]$to) == 1 && to %in% spec[[from]]$to) {
+              m[to, from] = paste("k", from, sep = "_")
+            }
+          }
+        }
+      }
+    } # }}}
+    model$coefmat <- m
+  }#}}}
+
+  # Try to create a function compiled from C code if more than one observed {{{
+  # variable and gcc is available
+  if (length(obs_vars) > 1) {
+    if (Sys.which("gcc") != "") {
+
+      # Translate the R code for the derivatives to C code
+      diffs.C <- paste(diffs, collapse = ";\n")
+      diffs.C <- paste0(diffs.C, ";")
+
+      # HS
+      diffs.C <- gsub(HS_decline, "(time <= tb ? k1 : k2)", diffs.C, fixed = TRUE)
+
+      for (i in seq_along(diffs)) {
+        state_var <- names(diffs)[i]
+
+        # IORE
+        if (state_var %in% obs_vars) {
+          if (spec[[state_var]]$type == "IORE") {
+            diffs.C <- gsub(paste0(state_var, "^N_", state_var),
+                            paste0("pow(y[", i - 1, "], N_", state_var, ")"),
+                            diffs.C, fixed = TRUE)
+          }
+        }
+
+        # Replace d_... terms by f[i-1]
+        # First line
+        pattern <- paste0("^d_", state_var)
+        replacement <- paste0("\nf[", i - 1, "]")
+        diffs.C <- gsub(pattern, replacement, diffs.C)
+        # Other lines
+        pattern <- paste0("\\nd_", state_var)
+        replacement <- paste0("\nf[", i - 1, "]")
+        diffs.C <- gsub(pattern, replacement, diffs.C)
+
+        # Replace names of observed variables by y[i],
+        # making the implicit assumption that the observed variables only occur after "* "
+        pattern <- paste0("\\* ", state_var)
+        replacement <- paste0("* y[", i - 1, "]")
+        diffs.C <- gsub(pattern, replacement, diffs.C)
+      }
+
+      derivs_sig <- signature(n = "integer", t = "numeric", y = "numeric",
+                              f = "numeric", rpar = "numeric", ipar = "integer")
+
+      # Declare the time variable in the body of the function if it is used
+      derivs_code <- if (spec[[1]]$type %in% c("FOMC", "DFOP", "HS")) {
+        paste0("double time = *t;\n", diffs.C)
+      } else {
+        diffs.C
+      }
+
+      # Define the function initializing the parameters
+      npar <- length(parms)
+      initpar_code <- paste0(
+        "static double parms [", npar, "];\n",
+        paste0("#define ", parms, " parms[", 0:(npar - 1), "]\n", collapse = ""),
+        "\n",
+        "void initpar(void (* odeparms)(int *, double *)) {\n",
+        "    int N = ", npar, ";\n",
+        "    odeparms(&N, parms);\n",
+        "}\n\n")
+
+      # Try to build a shared library
+      cf <- try(cfunction(list(func = derivs_sig), derivs_code,
+                          otherdefs = initpar_code,
+                          verbose = verbose,
+                          convention = ".C", language = "C"),
+                silent = TRUE)
+
+      if (!inherits(cf, "try-error")) {
+        if (!quiet) message("Successfully compiled differential equation model from auto-generated C code.")
+        model$cf <- cf
+      }
+    }
+  }
+  # }}}
+
+  class(model) <- "mkinmod"
+  return(model)
+}
+
+#' Print mkinmod objects
+#'
+#' Print mkinmod objects in a way that the user finds his way to get to its
+#' components.
+#'
+#' @param x An \code{\link{mkinmod}} object.
+#' @param \dots Not used.
+#' @examples
+#'
+#'   m_synth_SFO_lin <- mkinmod(parent = list(type = "SFO", to = "M1"),
+#'                              M1 = list(type = "SFO", to = "M2"),
+#'                              M2 = list(type = "SFO"), use_of_ff = "max")
+#'
+#'   print(m_synth_SFO_lin)
+#'
+#' @export
+print.mkinmod <- function(x, ...) {
+  cat("<mkinmod> model generated with\n")
+  cat("Use of formation fractions $use_of_ff:", x$use_of_ff, "\n")
+  cat("Specification $spec:\n")
+  for (obs in names(x$spec)) {
+    cat("$", obs, "\n", sep = "")
+    spl <- x$spec[[obs]]
+    cat("$type:", spl$type)
+    if (!is.null(spl$to) && length(spl$to)) cat("; $to: ", paste(spl$to, collapse = ", "), sep = "")
+    cat("; $sink: ", spl$sink, sep = "")
+    if (!is.null(spl$full_name)) if (!is.na(spl$full_name)) cat("; $full_name:", spl$full_name)
+    cat("\n")
+  }
+  if (is.matrix(x$coefmat)) cat("Coefficient matrix $coefmat available\n")
+  if (!is.null(x$cf)) cat("Compiled model $cf available\n")
+  cat("Differential equations:\n")
+  nice_diffs <- gsub("^(d.*) =", "\\1/dt =", x[["diffs"]])
+  writeLines(strwrap(nice_diffs, exdent = 11))
+}
+# vim: set foldmethod=marker ts=2 sw=2 expandtab: