jranke / mkin
1
#' Function to set up a kinetic model with one or more state variables
2
#'
3
#' This function is usually called using a call to [mkinsub()] for each observed
4
#' variable, specifying the corresponding submodel as well as outgoing pathways
5
#' (see examples).
6
#'
7
#' For the definition of model types and their parameters, the equations given
8
#' in the FOCUS and NAFTA guidance documents are used.
9
#'
10
#' For kinetic models with more than one observed variable, a symbolic solution
11
#' of the system of differential equations is included in the resulting
12
#' mkinmod object in some cases, speeding up the solution.
13
#'
14
#' If a C compiler is found by [pkgbuild::has_compiler()] and there
15
#' is more than one observed variable in the specification, C code is generated
16
#' for evaluating the differential equations, compiled using
17
#' [inline::cfunction()] and added to the resulting mkinmod object.
18
#'
19
#' @param ...  For each observed variable, a list as obtained by [mkinsub()]
20
#'   has to be specified as an argument (see examples).  Currently, single
21
#'   first order kinetics "SFO", indeterminate order rate equation kinetics
22
#'   "IORE", or single first order with reversible binding "SFORB" are
23
#'   implemented for all variables, while "FOMC", "DFOP", "HS" and "logistic"
24
#'   can additionally be chosen for the first variable which is assumed to be
25
#'   the source compartment.
26
#'   Additionally, [mkinsub()] has an argument \code{to}, specifying names of
27
#'   variables to which a transfer is to be assumed in the model.
28
#'   If the argument \code{use_of_ff} is set to "min"
29
#'   (default) and the model for the compartment is "SFO" or "SFORB", an
30
#'   additional [mkinsub()] argument can be \code{sink = FALSE}, effectively
31
#'   fixing the flux to sink to zero.
32
#'   In print.mkinmod, this argument is currently not used.
33
#' @param use_of_ff Specification of the use of formation fractions in the
34
#'   model equations and, if applicable, the coefficient matrix.  If "max",
35
#'   formation fractions are always used (default).  If "min", a minimum use of
36
#'   formation fractions is made, i.e. each first-order pathway to a metabolite
37
#'   has its own rate constant.
38
#' @param speclist The specification of the observed variables and their
39
#'   submodel types and pathways can be given as a single list using this
40
#'   argument. Default is NULL.
41
#' @param quiet Should messages be suppressed?
42
#' @param verbose If \code{TRUE}, passed to [inline::cfunction()] if
43
#'   applicable to give detailed information about the C function being built.
44
#' @param name A name for the model. Should be a valid R object name.
45
#' @param dll_dir Directory where an DLL object, if generated internally by
46
#'   [inline::cfunction()], should be saved.  The DLL will only be stored in a
47
#'   permanent location for use in future sessions, if 'dll_dir' and 'name'
48
#'   are specified.
49
#' @param unload If a DLL from the target location in 'dll_dir' is already
50
#'   loaded, should that be unloaded first?
51
#' @param overwrite If a file exists at the target DLL location in 'dll_dir',
52
#'   should this be overwritten?
53
#' @importFrom methods signature
54
#' @return A list of class \code{mkinmod} for use with [mkinfit()],
55
#'   containing, among others,
56
#'   \item{diffs}{
57
#'     A vector of string representations of differential equations, one for
58
#'     each modelling variable.
59
#'   }
60
#'   \item{map}{
61
#'     A list containing named character vectors for each observed variable,
62
#'     specifying the modelling variables by which it is represented.
63
#'   }
64
#'   \item{use_of_ff}{
65
#'     The content of \code{use_of_ff} is passed on in this list component.
66
#'   }
67
#'   \item{deg_func}{
68
#'     If generated, a function containing the solution of the degradation
69
#'     model.
70
#'   }
71
#'   \item{coefmat}{
72
#'     The coefficient matrix, if the system of differential equations can be
73
#'     represented by one.
74
#'   }
75
#'   \item{cf}{
76
#'     If generated, a compiled function calculating the derivatives as
77
#'     returned by cfunction.
78
#'   }
79
#' @note The IORE submodel is not well tested for metabolites. When using this
80
#'   model for metabolites, you may want to read the note in the help
81
#'   page to [mkinfit].
82
#' @author Johannes Ranke
83
#' @references FOCUS (2006) \dQuote{Guidance Document on Estimating Persistence
84
#'   and Degradation Kinetics from Environmental Fate Studies on Pesticides in
85
#'   EU Registration} Report of the FOCUS Work Group on Degradation Kinetics,
86
#'   EC Document Reference Sanco/10058/2005 version 2.0, 434 pp,
87
#'   \url{http://esdac.jrc.ec.europa.eu/projects/degradation-kinetics}
88
#'
89
#'   NAFTA Technical Working Group on Pesticides (not dated) Guidance for
90
#'   Evaluating and Calculating Degradation Kinetics in Environmental Media
91
#' @examples
92
#'
93
#' # Specify the SFO model (this is not needed any more, as we can now mkinfit("SFO", ...)
94
#' SFO <- mkinmod(parent = mkinsub("SFO"))
95
#'
96
#' # One parent compound, one metabolite, both single first order
97
#' SFO_SFO <- mkinmod(
98
#'   parent = mkinsub("SFO", "m1"),
99
#'   m1 = mkinsub("SFO"))
100
#' print(SFO_SFO)
101
#'
102
#' \dontrun{
103
#'  fit_sfo_sfo <- mkinfit(SFO_SFO, FOCUS_2006_D, quiet = TRUE, solution_type = "deSolve")
104
#'
105
#'  # Now supplying compound names used for plotting, and write to user defined location
106
#'  # We need to choose a path outside the session tempdir because this gets removed
107
#'  DLL_dir <- "~/.local/share/mkin"
108
#'  if (!dir.exists(DLL_dir)) dir.create(DLL_dir)
109
#'  SFO_SFO.2 <- mkinmod(
110
#'    parent = mkinsub("SFO", "m1", full_name = "Test compound"),
111
#'    m1 = mkinsub("SFO", full_name = "Metabolite M1"),
112
#'    name = "SFO_SFO", dll_dir = DLL_dir, unload = TRUE, overwrite = TRUE)
113
#' # Now we can save the model and restore it in a new session
114
#' saveRDS(SFO_SFO.2, file = "~/SFO_SFO.rds")
115
#' # Terminate the R session here if you would like to check, and then do
116
#' library(mkin)
117
#' SFO_SFO.3 <- readRDS("~/SFO_SFO.rds")
118
#' fit_sfo_sfo <- mkinfit(SFO_SFO.3, FOCUS_2006_D, quiet = TRUE, solution_type = "deSolve")
119
#'
120
#' # Show details of creating the C function
121
#' SFO_SFO <- mkinmod(
122
#'   parent = mkinsub("SFO", "m1"),
123
#'   m1 = mkinsub("SFO"), verbose = TRUE)
124
#'
125
#' # The symbolic solution which is available in this case is not
126
#' # made for human reading but for speed of computation
127
#' SFO_SFO$deg_func
128
#'
129
#' # If we have several parallel metabolites
130
#' # (compare tests/testthat/test_synthetic_data_for_UBA_2014.R)
131
#' m_synth_DFOP_par <- mkinmod(
132
#'  parent = mkinsub("DFOP", c("M1", "M2")),
133
#'  M1 = mkinsub("SFO"),
134
#'  M2 = mkinsub("SFO"),
135
#'  quiet = TRUE)
136
#'
137
#' fit_DFOP_par_c <- mkinfit(m_synth_DFOP_par,
138
#'   synthetic_data_for_UBA_2014[[12]]$data,
139
#'   quiet = TRUE)
140
#' }
141
#'
142
#' @export mkinmod
143
mkinmod <- function(..., use_of_ff = "max", name = NULL,
144
  speclist = NULL, quiet = FALSE, verbose = FALSE, dll_dir = NULL,
145
  unload = FALSE, overwrite = FALSE)
146
{
147 1
  if (is.null(speclist)) spec <- list(...)
148 1
  else spec <- speclist
149 1
  obs_vars <- names(spec)
150

151 1
  save_msg <- "You need to specify both 'name' and 'dll_dir' to save a model DLL"
152 1
  if (!is.null(dll_dir)) {
153 0
    if (!dir.exists(dll_dir)) stop(dll_dir, " does not exist")
154 0
    if (is.null(name)) stop(save_msg)
155
  }
156

157
  # Check if any of the names of the observed variables contains any other
158 1
  for (obs_var in obs_vars) {
159 1
    if (length(grep(obs_var, obs_vars)) > 1) stop("Sorry, variable names can not contain each other")
160 1
    if (grepl("_to_", obs_var)) stop("Sorry, names of observed variables can not contain _to_")
161 1
    if (obs_var == "sink") stop("Naming a compound 'sink' is not supported")
162
  }
163

164 1
  if (!use_of_ff %in% c("min", "max"))
165 1
    stop("The use of formation fractions 'use_of_ff' can only be 'min' or 'max'")
166

167 1
  parms <- vector()
168
  # }}}
169

170
  # Do not return a coefficient matrix mat when FOMC, IORE, DFOP, HS or logistic is used for the parent {{{
171 1
  if(spec[[1]]$type %in% c("FOMC", "IORE", "DFOP", "HS", "logistic")) {
172 1
    mat = FALSE
173 1
  } else mat = TRUE
174
  #}}}
175

176
  # Establish a list of differential equations as well as a map from observed {{{
177
  # compartments to differential equations
178 1
  diffs <- vector()
179 1
  map <- list()
180 1
  for (varname in obs_vars)
181
  {
182
    # Check the type component of the compartment specification {{{
183 0
    if(is.null(spec[[varname]]$type)) stop(
184 0
      "Every part of the model specification must be a list containing a type component")
185 1
    if(!spec[[varname]]$type %in% c("SFO", "FOMC", "IORE", "DFOP", "HS", "SFORB", "logistic")) stop(
186 1
      "Available types are SFO, FOMC, IORE, DFOP, HS, SFORB and logistic only")
187 1
    if(spec[[varname]]$type %in% c("FOMC", "DFOP", "HS", "logistic") & match(varname, obs_vars) != 1) {
188 1
        stop(paste("Types FOMC, DFOP, HS and logistic are only implemented for the first compartment,",
189 1
                   "which is assumed to be the source compartment"))
190
    }
191
    #}}}
192
    # New (sub)compartments (boxes) needed for the model type {{{
193 1
    new_boxes <- switch(spec[[varname]]$type,
194 1
      SFO = varname,
195 1
      FOMC = varname,
196 1
      IORE = varname,
197 1
      DFOP = varname,
198 1
      HS = varname,
199 1
      logistic = varname,
200 1
      SFORB = paste(varname, c("free", "bound"), sep = "_")
201
    )
202 1
    map[[varname]] <- new_boxes
203 1
    names(map[[varname]]) <- rep(spec[[varname]]$type, length(new_boxes)) #}}}
204
    # Start a new differential equation for each new box {{{
205 1
    new_diffs <- paste("d_", new_boxes, " =", sep = "")
206 1
    names(new_diffs) <- new_boxes
207 1
    diffs <- c(diffs, new_diffs) #}}}
208
  } #}}}
209

210
  # Create content of differential equations and build parameter list {{{
211 1
  for (varname in obs_vars)
212
  {
213
    # Get the name of the box(es) we are working on for the decline term(s)
214 1
    box_1 = map[[varname]][[1]] # This is the only box unless type is SFORB
215
    # Turn on sink if this is not explicitly excluded by the user by
216
    # specifying sink=FALSE
217 1
    if(is.null(spec[[varname]]$sink)) spec[[varname]]$sink <- TRUE
218 1
    if(spec[[varname]]$type %in% c("SFO", "IORE", "SFORB")) { # {{{ Add decline term
219 1
      if (use_of_ff == "min") { # Minimum use of formation fractions
220 1
        if(spec[[varname]]$type == "IORE" && length(spec[[varname]]$to) > 0) {
221 1
           stop("Transformation reactions from compounds modelled with IORE\n",
222 1
                "are only supported with formation fractions (use_of_ff = 'max')")
223
        }
224 1
        if(spec[[varname]]$sink) {
225
          # If sink is requested, add first-order/IORE sink term
226 1
          k_compound_sink <- paste("k", box_1, "sink", sep = "_")
227 1
          if(spec[[varname]]$type == "IORE") {
228 0
            k_compound_sink <- paste("k__iore", box_1, "sink", sep = "_")
229
          }
230 1
          parms <- c(parms, k_compound_sink)
231 1
          decline_term <- paste(k_compound_sink, "*", box_1)
232 1
          if(spec[[varname]]$type == "IORE") {
233 0
            N <- paste("N", box_1, sep = "_")
234 0
            parms <- c(parms, N)
235 0
            decline_term <- paste0(decline_term, "^", N)
236
          }
237
        } else { # otherwise no decline term needed here
238 1
          decline_term = "0"
239
        }
240
      } else { # Maximum use of formation fractions
241 1
        k_compound <- paste("k", box_1, sep = "_")
242 1
        if(spec[[varname]]$type == "IORE") {
243 1
          k_compound <- paste("k__iore", box_1, sep = "_")
244
        }
245 1
        parms <- c(parms, k_compound)
246 1
        decline_term <- paste(k_compound, "*", box_1)
247 1
        if(spec[[varname]]$type == "IORE") {
248 1
          N <- paste("N", box_1, sep = "_")
249 1
          parms <- c(parms, N)
250 1
          decline_term <- paste0(decline_term, "^", N)
251
        }
252
      }
253
    } #}}}
254 1
    if(spec[[varname]]$type == "FOMC") { # {{{ Add FOMC decline term
255
      # From p. 53 of the FOCUS kinetics report, without the power function so it works in C
256 1
      decline_term <- paste("(alpha/beta) * 1/((time/beta) + 1) *", box_1)
257 1
      parms <- c(parms, "alpha", "beta")
258
    } #}}}
259 1
    if(spec[[varname]]$type == "DFOP") { # {{{ Add DFOP decline term
260
      # From p. 57 of the FOCUS kinetics report
261 1
      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)
262 1
      parms <- c(parms, "k1", "k2", "g")
263
    } #}}}
264 1
    HS_decline <- "ifelse(time <= tb, k1, k2)" # Used below for automatic translation to C
265 1
    if(spec[[varname]]$type == "HS") { # {{{ Add HS decline term
266
      # From p. 55 of the FOCUS kinetics report
267 1
      decline_term <- paste(HS_decline, "*", box_1)
268 1
      parms <- c(parms, "k1", "k2", "tb")
269
    } #}}}
270 1
    if(spec[[varname]]$type == "logistic") { # {{{ Add logistic decline term
271
      # From p. 67 of the FOCUS kinetics report (2014)
272 1
      decline_term <- paste("(k0 * kmax)/(k0 + (kmax - k0) * exp(-r * time)) *", box_1)
273 1
      parms <- c(parms, "kmax", "k0", "r")
274
    } #}}}
275
    # Add origin decline term to box 1 (usually the only box, unless type is SFORB)#{{{
276 1
    diffs[[box_1]] <- paste(diffs[[box_1]], "-", decline_term)#}}}
277 1
    if(spec[[varname]]$type == "SFORB") { # {{{ Add SFORB reversible binding terms
278 1
      box_2 = map[[varname]][[2]]
279 1
      k_free_bound <- paste("k", varname, "free", "bound", sep = "_")
280 1
      k_bound_free <- paste("k", varname, "bound", "free", sep = "_")
281 1
      parms <- c(parms, k_free_bound, k_bound_free)
282 1
      reversible_binding_term_1 <- paste("-", k_free_bound, "*", box_1, "+",
283 1
        k_bound_free, "*", box_2)
284 1
      reversible_binding_term_2 <- paste("+", k_free_bound, "*", box_1, "-",
285 1
        k_bound_free, "*", box_2)
286 1
      diffs[[box_1]] <- paste(diffs[[box_1]], reversible_binding_term_1)
287 1
      diffs[[box_2]] <- paste(diffs[[box_2]], reversible_binding_term_2)
288
    } #}}}
289

290
    # Transfer between compartments#{{{
291 1
    to <- spec[[varname]]$to
292 1
    if(!is.null(to)) {
293
      # Name of box from which transfer takes place
294 1
      origin_box <- box_1
295

296
      # Number of targets
297 1
      n_targets = length(to)
298

299
      # Add transfer terms to listed compartments
300 1
      for (target in to) {
301 0
        if (!target %in% obs_vars) stop("You did not specify a submodel for target variable ", target)
302 1
        target_box <- switch(spec[[target]]$type,
303 1
          SFO = target,
304 1
          IORE = target,
305 1
          SFORB = paste(target, "free", sep = "_"))
306 1
        if (use_of_ff == "min" && spec[[varname]]$type %in% c("SFO", "SFORB"))
307
        {
308 1
          k_from_to <- paste("k", origin_box, target_box, sep = "_")
309 1
          parms <- c(parms, k_from_to)
310 1
          diffs[[origin_box]] <- paste(diffs[[origin_box]], "-",
311 1
            k_from_to, "*", origin_box)
312 1
          diffs[[target_box]] <- paste(diffs[[target_box]], "+",
313 1
            k_from_to, "*", origin_box)
314
        } else {
315
          # Do not introduce a formation fraction if this is the only target
316 1
          if (spec[[varname]]$sink == FALSE && n_targets == 1) {
317 1
            diffs[[target_box]] <- paste(diffs[[target_box]], "+",
318 1
                                         decline_term)
319
          } else {
320 1
            fraction_to_target = paste("f", origin_box, "to", target, sep = "_")
321 1
            parms <- c(parms, fraction_to_target)
322 1
            diffs[[target_box]] <- paste(diffs[[target_box]], "+",
323 1
                fraction_to_target, "*", decline_term)
324
          }
325
        }
326
      }
327
    } #}}}
328
  } #}}}
329

330 1
  model <- list(diffs = diffs, parms = parms, map = map, spec = spec, use_of_ff = use_of_ff, name = name)
331

332
  # Create coefficient matrix if possible #{{{
333 1
  if (mat) {
334 1
    boxes <- names(diffs)
335 1
    n <- length(boxes)
336 1
    m <- matrix(nrow=n, ncol=n, dimnames=list(boxes, boxes))
337

338 1
    if (use_of_ff == "min") { # {{{ Minimum use of formation fractions
339 1
      for (from in boxes) {
340 1
        for (to in boxes) {
341 1
          if (from == to) { # diagonal elements
342 1
            k.candidate = paste("k", from, c(boxes, "sink"), sep = "_")
343 1
            k.candidate = sub("free.*bound", "free_bound", k.candidate)
344 1
            k.candidate = sub("bound.*free", "bound_free", k.candidate)
345 1
            k.effective = intersect(model$parms, k.candidate)
346 1
            m[from,to] = ifelse(length(k.effective) > 0,
347 1
                paste("-", k.effective, collapse = " "), "0")
348

349
          } else {          # off-diagonal elements
350 1
            k.candidate = paste("k", from, to, sep = "_")
351 1
            if (sub("_free$", "", from) == sub("_bound$", "", to)) {
352 1
              k.candidate = paste("k", sub("_free$", "_free_bound", from), sep = "_")
353
            }
354 1
            if (sub("_bound$", "", from) == sub("_free$", "", to)) {
355 1
              k.candidate = paste("k", sub("_bound$", "_bound_free", from), sep = "_")
356
            }
357 1
            k.effective = intersect(model$parms, k.candidate)
358 1
            m[to, from] = ifelse(length(k.effective) > 0,
359 1
                k.effective, "0")
360
          }
361
        }
362
      } # }}}
363
    } else { # {{{ Use formation fractions where possible
364 1
      for (from in boxes) {
365 1
        for (to in boxes) {
366 1
          if (from == to) { # diagonal elements
367 1
            k.candidate = paste("k", from, sep = "_")
368 1
            m[from,to] = ifelse(k.candidate %in% model$parms,
369 1
                paste("-", k.candidate), "0")
370 1
            if(grepl("_free", from)) { # add transfer to bound compartment for SFORB
371 1
              m[from,to] = paste(m[from,to], "-", paste("k", from, "bound", sep = "_"))
372
            }
373 1
            if(grepl("_bound", from)) { # add backtransfer to free compartment for SFORB
374 1
              m[from,to] = paste("- k", from, "free", sep = "_")
375
            }
376 1
            m[from,to] = m[from,to]
377
          } else {          # off-diagonal elements
378 1
            f.candidate = paste("f", from, "to", to, sep = "_")
379 1
            k.candidate = paste("k", from, to, sep = "_")
380 1
            k.candidate = sub("free.*bound", "free_bound", k.candidate)
381 1
            k.candidate = sub("bound.*free", "bound_free", k.candidate)
382 1
            m[to, from] = ifelse(f.candidate %in% model$parms,
383 1
              paste(f.candidate, " * k_", from, sep = ""),
384 1
              ifelse(k.candidate %in% model$parms, k.candidate, "0"))
385
            # Special case: singular pathway and no sink
386 1
            if (spec[[from]]$sink == FALSE && length(spec[[from]]$to) == 1 && to %in% spec[[from]]$to) {
387 1
              m[to, from] = paste("k", from, sep = "_")
388
            }
389
          }
390
        }
391
      }
392
    } # }}}
393 1
    model$coefmat <- m
394
  }#}}}
395

396
  # Try to create a function compiled from C code if there is more than one observed variable {{{
397
  # and a compiler is available
398 1
  if (length(obs_vars) > 1 & pkgbuild::has_compiler()) {
399

400
    # Translate the R code for the derivatives to C code
401 1
    diffs.C <- paste(diffs, collapse = ";\n")
402 1
    diffs.C <- paste0(diffs.C, ";")
403

404
    # HS
405 1
    diffs.C <- gsub(HS_decline, "(time <= tb ? k1 : k2)", diffs.C, fixed = TRUE)
406

407 1
    for (i in seq_along(diffs)) {
408 1
      state_var <- names(diffs)[i]
409

410
      # IORE
411 1
      if (state_var %in% obs_vars) {
412 1
        if (spec[[state_var]]$type == "IORE") {
413 0
          diffs.C <- gsub(paste0(state_var, "^N_", state_var),
414 0
                          paste0("pow(y[", i - 1, "], N_", state_var, ")"),
415 0
                          diffs.C, fixed = TRUE)
416
        }
417
      }
418

419
      # Replace d_... terms by f[i-1]
420
      # First line
421 1
      pattern <- paste0("^d_", state_var)
422 1
      replacement <- paste0("\nf[", i - 1, "]")
423 1
      diffs.C <- gsub(pattern, replacement, diffs.C)
424
      # Other lines
425 1
      pattern <- paste0("\\nd_", state_var)
426 1
      replacement <- paste0("\nf[", i - 1, "]")
427 1
      diffs.C <- gsub(pattern, replacement, diffs.C)
428

429
      # Replace names of observed variables by y[i],
430
      # making the implicit assumption that the observed variables only occur after "* "
431 1
      pattern <- paste0("\\* ", state_var)
432 1
      replacement <- paste0("* y[", i - 1, "]")
433 1
      diffs.C <- gsub(pattern, replacement, diffs.C)
434
    }
435

436 1
    derivs_sig <- signature(n = "integer", t = "numeric", y = "numeric",
437 1
                            f = "numeric", rpar = "numeric", ipar = "integer")
438

439
    # Declare the time variable in the body of the function if it is used
440 1
    derivs_code <- if (spec[[1]]$type %in% c("FOMC", "DFOP", "HS")) {
441 1
      paste0("double time = *t;\n", diffs.C)
442
    } else {
443 1
      diffs.C
444
    }
445

446
    # Define the function initializing the parameters
447 1
    npar <- length(parms)
448 1
    initpar_code <- paste0(
449 1
      "static double parms [", npar, "];\n",
450 1
      paste0("#define ", parms, " parms[", 0:(npar - 1), "]\n", collapse = ""),
451 1
      "\n",
452 1
      "void initpar(void (* odeparms)(int *, double *)) {\n",
453 1
      "    int N = ", npar, ";\n",
454 1
      "    odeparms(&N, parms);\n",
455 1
      "}\n\n")
456

457
    # Try to build a shared library
458 1
    model$cf <- try(inline::cfunction(derivs_sig, derivs_code,
459 1
      otherdefs = initpar_code,
460 1
      verbose = verbose, name = "diffs",
461 1
      convention = ".C", language = "C"),
462 1
      silent = TRUE)
463

464 1
    if (!inherits(model$cf, "try-error")) {
465 1
      if (is.null(dll_dir)) {
466 0
        if (!quiet) message("Temporary DLL for differentials generated and loaded")
467 1
        dll_info <- inline::getDynLib(model$cf)
468
      } else {
469 0
        dll_info <- inline::moveDLL(model$cf, name, dll_dir,
470 0
          unload = unload, overwrite = overwrite, verbose = !quiet)
471
      }
472
    }
473
  }
474
  # }}}
475

476
  # Attach a degradation function if an analytical solution is available
477 1
  model$deg_func <- create_deg_func(spec, use_of_ff)
478

479 1
  class(model) <- "mkinmod"
480 1
  return(model)
481
}
482

483
#' Print mkinmod objects
484
#'
485
#' Print mkinmod objects in a way that the user finds his way to get to its
486
#' components.
487
#'
488
#' @rdname mkinmod
489
#' @param x An \code{\link{mkinmod}} object.
490
#' @export
491
print.mkinmod <- function(x, ...) {
492 1
  cat("<mkinmod> model generated with\n")
493 1
  cat("Use of formation fractions $use_of_ff:", x$use_of_ff, "\n")
494 1
  cat("Specification $spec:\n")
495 1
  for (obs in names(x$spec)) {
496 1
    cat("$", obs, "\n", sep = "")
497 1
    spl <- x$spec[[obs]]
498 1
    cat("$type:", spl$type)
499 1
    if (!is.null(spl$to) && length(spl$to)) cat("; $to: ", paste(spl$to, collapse = ", "), sep = "")
500 1
    cat("; $sink: ", spl$sink, sep = "")
501 0
    if (!is.null(spl$full_name)) if (!is.na(spl$full_name)) cat("; $full_name:", spl$full_name)
502 1
    cat("\n")
503
  }
504 1
  if (is.matrix(x$coefmat)) cat("Coefficient matrix $coefmat available\n")
505 0
  if (!is.null(x$cf)) cat("Compiled model $cf available\n")
506 1
  cat("Differential equations:\n")
507 1
  nice_diffs <- gsub("^(d.*) =", "\\1/dt =", x[["diffs"]])
508 1
  writeLines(strwrap(nice_diffs, exdent = 11))
509
}
510
# vim: set foldmethod=marker ts=2 sw=2 expandtab:

Read our documentation on viewing source code .

Loading