jranke / mkin
1
#' Fit one or more kinetic models with one or more state variables to one or
2
#' more datasets
3
#'
4
#' This function calls \code{\link{mkinfit}} on all combinations of models and
5
#' datasets specified in its first two arguments.
6
#'
7
#' @param models Either a character vector of shorthand names like
8
#'   \code{c("SFO", "FOMC", "DFOP", "HS", "SFORB")}, or an optionally named
9
#'   list of \code{\link{mkinmod}} objects.
10
#' @param datasets An optionally named list of datasets suitable as observed
11
#'   data for \code{\link{mkinfit}}.
12
#' @param cores The number of cores to be used for multicore processing. This
13
#'   is only used when the \code{cluster} argument is \code{NULL}. On Windows
14
#'   machines, cores > 1 is not supported, you need to use the \code{cluster}
15
#'   argument to use multiple logical processors. Per default, all cores
16
#'   detected by [parallel::detectCores()] are used, except on Windows where
17
#'   the default is 1.
18
#' @param cluster A cluster as returned by \code{\link{makeCluster}} to be used
19
#'   for parallel execution.
20
#' @param \dots Further arguments that will be passed to \code{\link{mkinfit}}.
21
#' @importFrom parallel mclapply parLapply detectCores
22
#' @return A two-dimensional \code{\link{array}} of \code{\link{mkinfit}}
23
#'   objects and/or try-errors that can be indexed using the model names for the
24
#'   first index (row index) and the dataset names for the second index (column
25
#'   index).
26
#' @author Johannes Ranke
27
#' @seealso \code{\link{[.mmkin}} for subsetting, \code{\link{plot.mmkin}} for
28
#'   plotting.
29
#' @keywords optimize
30
#' @examples
31
#'
32
#' \dontrun{
33
#' m_synth_SFO_lin <- mkinmod(parent = mkinsub("SFO", "M1"),
34
#'                            M1 = mkinsub("SFO", "M2"),
35
#'                            M2 = mkinsub("SFO"), use_of_ff = "max")
36
#'
37
#' m_synth_FOMC_lin <- mkinmod(parent = mkinsub("FOMC", "M1"),
38
#'                             M1 = mkinsub("SFO", "M2"),
39
#'                             M2 = mkinsub("SFO"), use_of_ff = "max")
40
#'
41
#' models <- list(SFO_lin = m_synth_SFO_lin, FOMC_lin = m_synth_FOMC_lin)
42
#' datasets <- lapply(synthetic_data_for_UBA_2014[1:3], function(x) x$data)
43
#' names(datasets) <- paste("Dataset", 1:3)
44
#'
45
#' time_default <- system.time(fits.0 <- mmkin(models, datasets, quiet = TRUE))
46
#' time_1 <- system.time(fits.4 <- mmkin(models, datasets, cores = 1, quiet = TRUE))
47
#'
48
#' time_default
49
#' time_1
50
#'
51
#' endpoints(fits.0[["SFO_lin", 2]])
52
#'
53
#' # plot.mkinfit handles rows or columns of mmkin result objects
54
#' plot(fits.0[1, ])
55
#' plot(fits.0[1, ], obs_var = c("M1", "M2"))
56
#' plot(fits.0[, 1])
57
#' # Use double brackets to extract a single mkinfit object, which will be plotted
58
#' # by plot.mkinfit and can be plotted using plot_sep
59
#' plot(fits.0[[1, 1]], sep_obs = TRUE, show_residuals = TRUE, show_errmin = TRUE)
60
#' plot_sep(fits.0[[1, 1]])
61
#' # Plotting with mmkin (single brackets, extracting an mmkin object) does not
62
#' # allow to plot the observed variables separately
63
#' plot(fits.0[1, 1])
64
#'
65
#' # On Windows, we can use multiple cores by making a cluster using the parallel
66
#' # package, which gets loaded with mkin, and passing it to mmkin, e.g.
67
#' cl <- makePSOCKcluster(12)
68
#' f <- mmkin(c("SFO", "FOMC", "DFOP"),
69
#'   list(A = FOCUS_2006_A, B = FOCUS_2006_B, C = FOCUS_2006_C, D = FOCUS_2006_D),
70
#'   cluster = cl, quiet = TRUE)
71
#' print(f)
72
#' # We get false convergence for the FOMC fit to FOCUS_2006_A because this
73
#' # dataset is really SFO, and the FOMC fit is overparameterised
74
#' stopCluster(cl)
75
#' }
76
#'
77
#' @export mmkin
78
mmkin <- function(models = c("SFO", "FOMC", "DFOP"), datasets,
79
  cores = if (Sys.info()["sysname"] == "Windows") 1 else parallel::detectCores(), cluster = NULL, ...)
80
{
81 1
  call <- match.call()
82 1
  parent_models_available = c("SFO", "FOMC", "DFOP", "HS", "SFORB", "IORE", "logistic")
83 1
  n.m <- length(models)
84 1
  n.d <- length(datasets)
85 1
  n.fits <- n.m * n.d
86 1
  fit_indices <- matrix(1:n.fits, ncol = n.d)
87

88
  # Check models and define their names
89 1
  if (!all(sapply(models, function(x) inherits(x, "mkinmod")))) {
90 1
    if (!all(models %in% parent_models_available)) {
91 0
      stop("Please supply models as a list of mkinmod objects or a vector combined of\n  ",
92 0
           paste(parent_models_available, collapse = ", "))
93
    } else {
94 1
      names(models) <- models
95
    }
96
  } else {
97 0
    if (is.null(names(models))) names(models) <- as.character(1:n.m)
98
  }
99

100
  # Check datasets and define their names
101 1
  if (is.null(names(datasets))) names(datasets) <- as.character(1:n.d)
102

103
  # Define names for fit index
104 1
  dimnames(fit_indices) <- list(model = names(models),
105 1
                                dataset = names(datasets))
106

107

108 1
  fit_function <- function(fit_index) {
109 1
    w <- which(fit_indices == fit_index, arr.ind = TRUE)
110 1
    model_index <- w[1]
111 1
    dataset_index <- w[2]
112 1
    res <- try(mkinfit(models[[model_index]], datasets[[dataset_index]], ...))
113 1
    if (!inherits(res, "try-error")) res$mkinmod$name <- names(models)[model_index]
114 1
    return(res)
115
  }
116

117 1
  if (is.null(cluster)) {
118 1
    results <- parallel::mclapply(as.list(1:n.fits), fit_function,
119 1
      mc.cores = cores, mc.preschedule = FALSE)
120
  } else {
121 0
    results <- parallel::parLapply(cluster, as.list(1:n.fits), fit_function)
122
  }
123

124 1
  attributes(results) <- attributes(fit_indices)
125 1
  attr(results, "call") <- call
126 1
  class(results) <- "mmkin"
127 1
  return(results)
128
}
129

130
#' Subsetting method for mmkin objects
131
#'
132
#' @param x An \code{\link{mmkin} object}
133
#' @param i Row index selecting the fits for specific models
134
#' @param j Column index selecting the fits to specific datasets
135
#' @param ... Not used, only there to satisfy the generic method definition
136
#' @param drop If FALSE, the method always returns an mmkin object, otherwise
137
#'   either a list of mkinfit objects or a single mkinfit object.
138
#' @return An object of class \code{\link{mmkin}}.
139
#' @author Johannes Ranke
140
#' @rdname Extract.mmkin
141
#' @examples
142
#'
143
#'   # Only use one core, to pass R CMD check --as-cran
144
#'   fits <- mmkin(c("SFO", "FOMC"), list(B = FOCUS_2006_B, C = FOCUS_2006_C),
145
#'                 cores = 1, quiet = TRUE)
146
#'   fits["FOMC", ]
147
#'   fits[, "B"]
148
#'   fits["SFO", "B"]
149
#'
150
#'   head(
151
#'     # This extracts an mkinfit object with lots of components
152
#'     fits[["FOMC", "B"]]
153
#'   )
154
#' @export
155
`[.mmkin` <- function(x, i, j, ..., drop = FALSE) {
156 1
  class(x) <- NULL
157 1
  x_sub <- x[i, j, drop = drop]
158 1
  if (!drop) class(x_sub) <- "mmkin"
159 1
  return(x_sub)
160
}
161

162
#' Print method for mmkin objects
163
#'
164
#' @param x An [mmkin] object.
165
#' @param \dots Not used.
166
#' @rdname mmkin
167
#' @export
168
print.mmkin <- function(x, ...) {
169 1
  cat("<mmkin> object\n")
170 1
  cat("Status of individual fits:\n\n")
171 1
  all_summary_warnings <- character()
172 1
  sww <- 0 # Counter for Shapiro-Wilks warnings
173

174 1
  display <- lapply(x,
175 1
    function(fit) {
176 0
      if (inherits(fit, "try-error")) return("E")
177 1
      sw <- fit$summary_warnings
178 1
      swn <- names(sw)
179 1
      if (length(sw) > 0) {
180 0
        if (any(grepl("S", swn))) {
181 0
          sww <<- sww + 1
182 0
          swn <- gsub("S", paste0("S", sww), swn)
183
        }
184 0
        warnstring <- paste(swn, collapse = ", ")
185 0
        names(sw) <- swn
186 0
        all_summary_warnings <<- c(all_summary_warnings, sw)
187 0
        return(warnstring)
188
      } else {
189 1
        return("OK")
190
      }
191
    })
192 1
  display <- unlist(display)
193 1
  dim(display) <- dim(x)
194 1
  dimnames(display) <- dimnames(x)
195 1
  print(display, quote = FALSE)
196

197 1
  cat("\n")
198 1
  if (any(display == "OK")) cat("OK: No warnings\n")
199 0
  if (any(display == "E")) cat("E: Error\n")
200 1
  u_swn <- unique(names(all_summary_warnings))
201 1
  u_w <- all_summary_warnings[u_swn]
202 1
  for (i in seq_along(u_w)) {
203 0
    cat(names(u_w)[i], ": ", u_w[i], "\n", sep = "")
204
  }
205

206
}
207

208
#' @export
209
update.mmkin <- function(object, ..., evaluate = TRUE)
210
{
211 0
  call <- attr(object, "call")
212

213 0
  update_arguments <- match.call(expand.dots = FALSE)$...
214

215 0
  if (length(update_arguments) > 0) {
216 0
    update_arguments_in_call <- !is.na(match(names(update_arguments), names(call)))
217
  }
218

219 0
  for (a in names(update_arguments)[update_arguments_in_call]) {
220 0
    call[[a]] <- update_arguments[[a]]
221
  }
222

223 0
  update_arguments_not_in_call <- !update_arguments_in_call
224 0
  if(any(update_arguments_not_in_call)) {
225 0
    call <- c(as.list(call), update_arguments[update_arguments_not_in_call])
226 0
    call <- as.call(call)
227
  }
228

229 0
  if(evaluate) eval(call, parent.frame())
230 0
  else call
231
}

Read our documentation on viewing source code .

Loading