From f3f415520c89f9d8526bf6fadc862ebd44be220d Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Thu, 17 Nov 2016 18:14:32 +0100 Subject: Remove trailing whitespace, clean headers Also ignore test.R in the top level directory, as it is not meant to be public --- R/SFORB.solution.R | 4 +-- R/add_err.R | 6 ++-- R/endpoints.R | 12 +++---- R/ilr.R | 4 +-- R/mkin_long_to_wide.R | 4 +-- R/mkinds.R | 2 +- R/mkinerrmin.R | 14 ++++---- R/mkinfit.R | 90 +++++++++++++++++++++++++------------------------- R/mkinmod.R | 20 +++++------ R/mkinparplot.R | 22 ++++++------ R/mkinpredict.R | 26 +++++++-------- R/mkinresplot.R | 10 +++--- R/mmkin.R | 10 +++--- R/plot.mkinfit.R | 32 +++++++++--------- R/plot.mmkin.R | 10 +++--- R/transform_odeparms.R | 10 +++--- 16 files changed, 136 insertions(+), 140 deletions(-) (limited to 'R') diff --git a/R/SFORB.solution.R b/R/SFORB.solution.R index 45a4533a..4cb94def 100644 --- a/R/SFORB.solution.R +++ b/R/SFORB.solution.R @@ -2,8 +2,8 @@ SFORB.solution = function(t, parent.0, k_12, k_21, k_1output) { sqrt_exp = sqrt(1/4 * (k_12 + k_21 + k_1output)^2 + k_12 * k_21 - (k_12 + k_1output) * k_21) b1 = 0.5 * (k_12 + k_21 + k_1output) + sqrt_exp b2 = 0.5 * (k_12 + k_21 + k_1output) - sqrt_exp - - parent = parent.0 * + + parent = parent.0 * (((k_12 + k_21 - b1)/(b2 - b1)) * exp(-b1 * t) + ((k_12 + k_21 - b2)/(b1 - b2)) * exp(-b2 * t)) } diff --git a/R/add_err.R b/R/add_err.R index 0995e634..4d998e94 100644 --- a/R/add_err.R +++ b/R/add_err.R @@ -16,8 +16,8 @@ # You should have received a copy of the GNU General Public License along with # this program. If not, see -add_err = function(prediction, sdfunc, - n = 1000, LOD = 0.1, reps = 2, +add_err = function(prediction, sdfunc, + n = 1000, LOD = 0.1, reps = 2, digits = 1, seed = NA) { if (!is.na(seed)) set.seed(seed) @@ -32,7 +32,7 @@ add_err = function(prediction, sdfunc, for (i in 1:n) { d_rep = data.frame(lapply(d_long, rep, each = 2)) d_rep$value = rnorm(length(d_rep$value), d_rep$value, sdfunc(d_rep$value)) - + d_rep[d_rep$time == 0 & d_rep$name %in% c("M1", "M2"), "value"] <- 0 # Set values below the LOD to NA diff --git a/R/endpoints.R b/R/endpoints.R index 03c38ee0..ac1e3e7c 100644 --- a/R/endpoints.R +++ b/R/endpoints.R @@ -8,19 +8,19 @@ endpoints <- function(fit) { parms.all <- c(fit$bparms.optim, fit$bparms.fixed) ep$ff <- vector() ep$SFORB <- vector() - ep$distimes <- data.frame(DT50 = rep(NA, length(obs_vars)), - DT90 = rep(NA, length(obs_vars)), + ep$distimes <- data.frame(DT50 = rep(NA, length(obs_vars)), + DT90 = rep(NA, length(obs_vars)), row.names = obs_vars) for (obs_var in obs_vars) { - type = names(fit$mkinmod$map[[obs_var]])[1] + type = names(fit$mkinmod$map[[obs_var]])[1] # Get formation fractions if directly fitted, and calculate remaining fraction to sink f_names = grep(paste("^f", obs_var, sep = "_"), names(parms.all), value=TRUE) if (length(f_names) > 0) { f_values = parms.all[f_names] f_to_sink = 1 - sum(f_values) - names(f_to_sink) = ifelse(type == "SFORB", - paste(obs_var, "free", "sink", sep = "_"), + names(f_to_sink) = ifelse(type == "SFORB", + paste(obs_var, "free", "sink", sep = "_"), paste(obs_var, "sink", sep = "_")) for (f_name in f_names) { ep$ff[[sub("f_", "", sub("_to_", "_", f_name))]] = f_values[[f_name]] @@ -91,7 +91,7 @@ endpoints <- function(fit) { silent = TRUE) if (inherits(DT50, "try-error")) DT50 = NA if (inherits(DT90, "try-error")) DT90 = NA - + ep$distimes[obs_var, c("DT50_k1")] = DT50_k1 ep$distimes[obs_var, c("DT50_k2")] = DT50_k2 } diff --git a/R/ilr.R b/R/ilr.R index 389653e7..620afc49 100644 --- a/R/ilr.R +++ b/R/ilr.R @@ -1,5 +1,3 @@ -# $Id$ - # Copyright (C) 2012 René Lehmann, Johannes Ranke # Contact: jranke@uni-bremen.de @@ -44,7 +42,7 @@ invilr<-function(x) { # Work around a numerical problem with NaN values returned by the above # Only works if there is only one NaN value: replace it with 1 # if the sum of the other components is < 1e-10 - if (sum(is.na(z)) == 1 && sum(z[!is.na(z)]) < 1e-10) + if (sum(is.na(z)) == 1 && sum(z[!is.na(z)]) < 1e-10) z = ifelse(is.na(z), 1, z) return(z) diff --git a/R/mkin_long_to_wide.R b/R/mkin_long_to_wide.R index 87e1afd5..081262f8 100644 --- a/R/mkin_long_to_wide.R +++ b/R/mkin_long_to_wide.R @@ -1,7 +1,5 @@ -# $Id$ - # Copyright (C) 2010-2011 Johannes Ranke -# Contact: mkin-devel@lists.berlios.de +# Contact: jranke@uni-bremen.de # This file is part of the R package mkin diff --git a/R/mkinds.R b/R/mkinds.R index cb3c712a..60924c21 100644 --- a/R/mkinds.R +++ b/R/mkinds.R @@ -31,7 +31,7 @@ #' @field data A dataframe with at least the columns name, time and value #' in order to be compatible with mkinfit #' @example inst/examples/mkinds.R -mkinds <- R6Class("mkinds", +mkinds <- R6Class("mkinds", public = list( title = NULL, sampling_times = NULL, diff --git a/R/mkinerrmin.R b/R/mkinerrmin.R index 6103dd1d..b361c466 100644 --- a/R/mkinerrmin.R +++ b/R/mkinerrmin.R @@ -24,7 +24,7 @@ mkinerrmin <- function(fit, alpha = 0.05) kinerrmin <- function(errdata, n.parms) { means.mean <- mean(errdata$value_mean, na.rm = TRUE) df = length(errdata$value_mean) - n.parms - + err.min <- sqrt((1 / qchisq(1 - alpha, df)) * sum((errdata$value_mean - errdata$value_pred)^2)/(means.mean^2)) @@ -32,12 +32,12 @@ mkinerrmin <- function(fit, alpha = 0.05) } means <- aggregate(value ~ time + name, data = fit$observed, mean, na.rm=TRUE) - errdata <- merge(means, fit$predicted, by = c("time", "name"), + errdata <- merge(means, fit$predicted, by = c("time", "name"), suffixes = c("_mean", "_pred")) errdata <- errdata[order(errdata$time, errdata$name), ] # Remove values at time zero for variables whose value for state.ini is fixed, - # as these will not have any effect in the optimization and should therefore not + # as these will not have any effect in the optimization and should therefore not # be counted as degrees of freedom. fixed_initials = gsub("_0$", "", rownames(subset(fit$fixed, type = "state"))) errdata <- subset(errdata, !(time == 0 & name %in% fixed_initials)) @@ -59,10 +59,10 @@ mkinerrmin <- function(fit, alpha = 0.05) # Rate constants and IORE exponents are attributed to the source variable n.k.optim <- length(grep(paste("^k", obs_var, sep="_"), names(parms.optim))) - n.k.optim <- n.k.optim + length(grep(paste("^log_k", obs_var, sep="_"), + n.k.optim <- n.k.optim + length(grep(paste("^log_k", obs_var, sep="_"), names(parms.optim))) n.k__iore.optim <- length(grep(paste("^k__iore", obs_var, sep="_"), names(parms.optim))) - n.k__iore.optim <- n.k__iore.optim + length(grep(paste("^log_k__iore", obs_var, + n.k__iore.optim <- n.k__iore.optim + length(grep(paste("^log_k__iore", obs_var, sep = "_"), names(parms.optim))) @@ -87,8 +87,8 @@ mkinerrmin <- function(fit, alpha = 0.05) # FOMC, DFOP and HS parameters are only counted if we are looking at the # first variable in the model which is always the source variable if (obs_var == fit$obs_vars[[1]]) { - special_parms = c("alpha", "log_alpha", "beta", "log_beta", - "k1", "log_k1", "k2", "log_k2", + special_parms = c("alpha", "log_alpha", "beta", "log_beta", + "k1", "log_k1", "k2", "log_k2", "g", "g_ilr", "tb", "log_tb") n.optim <- n.optim + length(intersect(special_parms, names(parms.optim))) } diff --git a/R/mkinfit.R b/R/mkinfit.R index 4c4a9717..05986ff9 100644 --- a/R/mkinfit.R +++ b/R/mkinfit.R @@ -23,7 +23,7 @@ if(getRversion() >= '2.15.1') utils::globalVariables(c("name", "time", "value")) mkinfit <- function(mkinmod, observed, parms.ini = "auto", - state.ini = "auto", + state.ini = "auto", fixed_parms = NULL, fixed_initials = names(mkinmod$diffs)[-1], from_max_mean = FALSE, @@ -45,7 +45,7 @@ mkinfit <- function(mkinmod, observed, { # Check mkinmod and generate a model for the variable whith the highest value # if a suitable string is given - parent_models_available = c("SFO", "FOMC", "DFOP", "HS", "SFORB", "IORE") + parent_models_available = c("SFO", "FOMC", "DFOP", "HS", "SFORB", "IORE") if (class(mkinmod) != "mkinmod") { presumed_parent_name = observed[which.max(observed$value), "name"] if (mkinmod[[1]] %in% parent_models_available) { @@ -55,7 +55,7 @@ mkinfit <- function(mkinmod, observed, } else { stop("Argument mkinmod must be of class mkinmod or a string containing one of\n ", paste(parent_models_available, collapse = ", ")) - } + } } # Check optimisation method and set maximum number of iterations if specified @@ -83,7 +83,7 @@ mkinfit <- function(mkinmod, observed, # Obtain data for decline from maximum mean value if requested if (from_max_mean) { # This is only used for simple decline models - if (length(obs_vars) > 1) + if (length(obs_vars) > 1) stop("Decline from maximum is only implemented for models with a single observed variable") means <- aggregate(value ~ time, data = observed, mean, na.rm=TRUE) @@ -138,7 +138,7 @@ mkinfit <- function(mkinmod, observed, k_salt = k_salt + 1e-4 } # Default values for rate constants for reversible binding - if (grepl("free_bound$", parmname)) parms.ini[parmname] = 0.1 + if (grepl("free_bound$", parmname)) parms.ini[parmname] = 0.1 if (grepl("bound_free$", parmname)) parms.ini[parmname] = 0.02 # Default values for IORE exponents if (grepl("^N", parmname)) parms.ini[parmname] = 1 @@ -282,9 +282,9 @@ mkinfit <- function(mkinmod, observed, transform_fractions = transform_fractions) # Solve the system with current transformed parameter values - out <- mkinpredict(mkinmod, parms, - odeini, outtimes, - solution_type = solution_type, + out <- mkinpredict(mkinmod, parms, + odeini, outtimes, + solution_type = solution_type, use_compiled = use_compiled, method.ode = method.ode, atol = atol, rtol = rtol, ...) @@ -303,27 +303,27 @@ mkinfit <- function(mkinmod, observed, outtimes_plot = seq(min(observed$time), max(observed$time), length.out=100) out_plot <- mkinpredict(mkinmod, parms, - odeini, outtimes_plot, - solution_type = solution_type, + odeini, outtimes_plot, + solution_type = solution_type, use_compiled = use_compiled, method.ode = method.ode, atol = atol, rtol = rtol, ...) - plot(0, type="n", + plot(0, type="n", xlim = range(observed$time), ylim = c(0, max(observed$value, na.rm=TRUE)), xlab = "Time", ylab = "Observed") col_obs <- pch_obs <- 1:length(obs_vars) lty_obs <- rep(1, length(obs_vars)) names(col_obs) <- names(pch_obs) <- names(lty_obs) <- obs_vars for (obs_var in obs_vars) { - points(subset(observed, name == obs_var, c(time, value)), + points(subset(observed, name == obs_var, c(time, value)), pch = pch_obs[obs_var], col = col_obs[obs_var]) } matlines(out_plot$time, out_plot[-1], col = col_obs, lty = lty_obs) - legend("topright", inset=c(0.05, 0.05), legend=obs_vars, + legend("topright", inset=c(0.05, 0.05), legend=obs_vars, col=col_obs, pch=pch_obs, lty=1:length(pch_obs)) } - + assign("cost.old", mC$model, inherits=TRUE) } return(mC) @@ -343,9 +343,9 @@ mkinfit <- function(mkinmod, observed, odeparms <- c(P[(length(state.ini.optim) + 1):length(P)], parms.fixed) # Solve the system with current transformed parameter values - out <- mkinpredict(mkinmod, odeparms, - odeini, outtimes, - solution_type = solution_type, + out <- mkinpredict(mkinmod, odeparms, + odeini, outtimes, + solution_type = solution_type, use_compiled = use_compiled, method.ode = method.ode, atol = atol, rtol = rtol, ...) @@ -381,8 +381,8 @@ mkinfit <- function(mkinmod, observed, # Do the fit and take the time fit_time <- system.time({ - fit <- modFit(cost, c(state.ini.optim, transparms.optim), - method = method.modFit, control = control.modFit, + fit <- modFit(cost, c(state.ini.optim, transparms.optim), + method = method.modFit, control = control.modFit, lower = lower, upper = upper, ...) # Reiterate the fit until convergence of the variance components (IRLS) @@ -423,7 +423,7 @@ mkinfit <- function(mkinmod, observed, # Check for convergence if (method.modFit == "Marq") { if (!fit$info %in% c(1, 2, 3)) { - fit$warning = paste0("Optimisation by method ", method.modFit, + fit$warning = paste0("Optimisation by method ", method.modFit, " did not converge.\n", "The message returned by nls.lm is:\n", fit$message) @@ -435,10 +435,10 @@ mkinfit <- function(mkinmod, observed, } if (method.modFit %in% c("Port", "SANN", "Nelder-Mead", "BFGS", "CG", "L-BFGS-B")) { if (fit$convergence != 0) { - fit$warning = paste0("Optimisation by method ", method.modFit, + fit$warning = paste0("Optimisation by method ", method.modFit, " did not converge.\n", "Convergence code is ", fit$convergence, - ifelse(is.null(fit$message), "", + ifelse(is.null(fit$message), "", paste0("\nMessage is ", fit$message))) warning(fit$warning) } else { @@ -483,7 +483,7 @@ mkinfit <- function(mkinmod, observed, # Estimate the Hessian for the model cost without parameter transformations # to make it possible to obtain the usual t-test # Code ported from FME::modFit - Jac_notrans <- try(gradient(function(p, ...) cost_notrans(p)$residuals$res, + Jac_notrans <- try(gradient(function(p, ...) cost_notrans(p)$residuals$res, bparms.optim, centered = TRUE), silent = TRUE) if (inherits(Jac_notrans, "try-error")) { warning("Calculation of the Jacobian failed for the cost function of the untransformed model.\n", @@ -494,9 +494,9 @@ mkinfit <- function(mkinmod, observed, } # Collect initial parameter values in three dataframes - fit$start <- data.frame(value = c(state.ini.optim, + fit$start <- data.frame(value = c(state.ini.optim, parms.optim)) - fit$start$type = c(rep("state", length(state.ini.optim)), + fit$start$type = c(rep("state", length(state.ini.optim)), rep("deparm", length(parms.optim))) fit$start_transformed = data.frame( @@ -505,7 +505,7 @@ mkinfit <- function(mkinmod, observed, upper = upper) fit$fixed <- data.frame(value = c(state.ini.fixed, parms.fixed)) - fit$fixed$type = c(rep("state", length(state.ini.fixed)), + fit$fixed$type = c(rep("state", length(state.ini.fixed)), rep("deparm", length(parms.fixed))) # Collect observed, predicted, residuals and weighting @@ -535,18 +535,18 @@ mkinfit <- function(mkinmod, observed, fit$reweight.max.iter <- reweight.max.iter # Return different sets of backtransformed parameters for summary and plotting - fit$bparms.optim <- bparms.optim + fit$bparms.optim <- bparms.optim fit$bparms.fixed <- bparms.fixed # Return ode and state parameters for further fitting - fit$bparms.ode <- bparms.all[mkinmod$parms] + fit$bparms.ode <- bparms.all[mkinmod$parms] fit$bparms.state <- c(bparms.all[setdiff(names(bparms.all), names(fit$bparms.ode))], state.ini.fixed) names(fit$bparms.state) <- gsub("_0$", "", names(fit$bparms.state)) fit$date <- date() - class(fit) <- c("mkinfit", "modFit") + class(fit) <- c("mkinfit", "modFit") return(fit) } @@ -602,36 +602,36 @@ summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, alpha = 0.05, par.lower <- param[pname, "Lower"] par.upper <- param[pname, "Upper"] names(par.lower) <- names(par.upper) <- pname - bpl <- backtransform_odeparms(par.lower, object$mkinmod, - object$transform_rates, + bpl <- backtransform_odeparms(par.lower, object$mkinmod, + object$transform_rates, object$transform_fractions) bpu <- backtransform_odeparms(par.upper, object$mkinmod, - object$transform_rates, + object$transform_rates, object$transform_fractions) bparam[names(bpl), "Lower"] <- bpl bparam[names(bpu), "Upper"] <- bpu - } + } } ans <- list( version = as.character(utils::packageVersion("mkin")), Rversion = paste(R.version$major, R.version$minor, sep="."), - date.fit = object$date, - date.summary = date(), - solution_type = object$solution_type, - method.modFit = object$method.modFit, - warning = object$warning, - use_of_ff = object$mkinmod$use_of_ff, + date.fit = object$date, + date.summary = date(), + solution_type = object$solution_type, + method.modFit = object$method.modFit, + warning = object$warning, + use_of_ff = object$mkinmod$use_of_ff, weight.ini = object$weight.ini, reweight.method = object$reweight.method, residuals = object$residuals, residualVariance = resvar, sigma = sqrt(resvar), modVariance = modVariance, - df = c(p, rdf), + df = c(p, rdf), cov.unscaled = covar, cov.scaled = covar * resvar, - info = object$info, + info = object$info, niter = object$iterations, calls = object$calls, time = object$time, @@ -654,8 +654,8 @@ summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, alpha = 0.05, ans$ff <- ep$ff if(distimes) ans$distimes <- ep$distimes if(length(ep$SFORB) != 0) ans$SFORB <- ep$SFORB - class(ans) <- c("summary.mkinfit", "summary.modFit") - return(ans) + class(ans) <- c("summary.mkinfit", "summary.modFit") + return(ans) } # Expanded from print.summary.modFit @@ -674,7 +674,7 @@ print.summary.mkinfit <- function(x, digits = max(3, getOption("digits") - 3), . cat("\nModel predictions using solution type", x$solution_type, "\n") - cat("\nFitted with method", x$method.modFit, + cat("\nFitted with method", x$method.modFit, "using", x$calls, "model solutions performed in", x$time[["elapsed"]], "s\n") cat("\nWeighting:", x$weight.ini) @@ -691,7 +691,7 @@ print.summary.mkinfit <- function(x, digits = max(3, getOption("digits") - 3), . cat("\nFixed parameter values:\n") if(length(x$fixed$value) == 0) cat("None\n") else print(x$fixed) - + cat("\nOptimised, transformed parameters with symmetric confidence intervals:\n") print(signif(x$par, digits = digits)) diff --git a/R/mkinmod.R b/R/mkinmod.R index 5ec62bbe..7f7f8587 100644 --- a/R/mkinmod.R +++ b/R/mkinmod.R @@ -44,7 +44,7 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL, quiet = FALSE, verb # Do not return a coefficient matrix mat when FOMC, IORE, DFOP or HS is used for the parent {{{ if(spec[[1]]$type %in% c("FOMC", "IORE", "DFOP", "HS")) { - mat = FALSE + mat = FALSE } else mat = TRUE #}}} @@ -60,7 +60,7 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL, quiet = FALSE, verb if(!spec[[varname]]$type %in% c("SFO", "FOMC", "IORE", "DFOP", "HS", "SFORB")) stop( "Available types are SFO, FOMC, IORE, 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,", + stop(paste("Types FOMC, DFOP and HS are only implemented for the first compartment,", "which is assumed to be the source compartment")) } #}}} @@ -109,7 +109,7 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL, quiet = FALSE, verb decline_term <- paste0(decline_term, "^", N) } } else { # otherwise no decline term needed here - decline_term = "0" + decline_term = "0" } } else { k_compound <- paste("k", box_1, sep = "_") @@ -156,7 +156,7 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL, quiet = FALSE, verb } 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 + # 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 = "_") @@ -189,9 +189,9 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL, quiet = FALSE, verb { k_from_to <- paste("k", origin_box, target_box, sep = "_") parms <- c(parms, k_from_to) - diffs[[origin_box]] <- paste(diffs[[origin_box]], "-", + diffs[[origin_box]] <- paste(diffs[[origin_box]], "-", k_from_to, "*", origin_box) - diffs[[target_box]] <- paste(diffs[[target_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 @@ -201,7 +201,7 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL, quiet = FALSE, verb } else { fraction_to_target = paste("f", origin_box, "to", target, sep = "_") parms <- c(parms, fraction_to_target) - diffs[[target_box]] <- paste(diffs[[target_box]], "+", + diffs[[target_box]] <- paste(diffs[[target_box]], "+", fraction_to_target, "*", decline_term) } } @@ -261,7 +261,7 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL, quiet = FALSE, verb 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 = ""), + 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) { @@ -317,7 +317,7 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL, quiet = FALSE, verb 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) @@ -369,6 +369,6 @@ print.mkinmod <- function(x, ...) { cat("\n") } if (is.matrix(x$coefmat)) cat("Coefficient matrix $coefmat available\n") - if (!is.null(x$cf)) cat("Compiled model $cf available\n") + if (!is.null(x$cf)) cat("Compiled model $cf available\n") } # vim: set foldmethod=marker ts=2 sw=2 expandtab: diff --git a/R/mkinparplot.R b/R/mkinparplot.R index 0e8bdf5a..af28e3a8 100644 --- a/R/mkinparplot.R +++ b/R/mkinparplot.R @@ -23,8 +23,8 @@ mkinparplot <- function(object) { if ("g" %in% deparms.optim) fractions.optim <- c("g", fractions.optim) rates.optim.unsorted = setdiff(deparms.optim, union(fractions.optim, N.optim)) rates.optim <- rownames(object$start[rates.optim.unsorted, ]) - n.plot <- c(state.optim = length(state.optim), - rates.optim = length(rates.optim), + n.plot <- c(state.optim = length(state.optim), + rates.optim = length(rates.optim), N.optim = length(N.optim), fractions.optim = length(fractions.optim)) n.plot <- n.plot[n.plot > 0] @@ -41,31 +41,31 @@ mkinparplot <- function(object) { values <- bpar[parnames] values_with_confints <- data.frame(t(subset(data.frame(t(values)), !is.na("Lower")))) xlim = switch(type, - state.optim = range(c(0, unlist(values)), + state.optim = range(c(0, unlist(values)), na.rm = TRUE, finite = TRUE), - rates.optim = range(c(0, unlist(values)), + rates.optim = range(c(0, unlist(values)), na.rm = TRUE, finite = TRUE), - N.optim = range(c(0, 1, unlist(values)), + N.optim = range(c(0, 1, unlist(values)), na.rm = TRUE, finite = TRUE), - fractions.optim = range(c(0, 1, unlist(values)), + fractions.optim = range(c(0, 1, unlist(values)), na.rm = TRUE, finite = TRUE)) parname_index <- length(parnames):1 # Reverse order for strip chart - stripchart(values["Estimate", ][parname_index], + stripchart(values["Estimate", ][parname_index], xlim = xlim, ylim = c(0.5, length(get(type)) + 0.5), yaxt = "n") if (type %in% c("rates.optim", "fractions.optim")) abline(v = 0, lty = 2) if (type %in% c("N.optim", "fractions.optim")) abline(v = 1, lty = 2) position <- ifelse(values["Estimate", ] < mean(xlim), "right", "left") - text(ifelse(position == "left", min(xlim), max(xlim)), - parname_index, parnames, + text(ifelse(position == "left", min(xlim), max(xlim)), + parname_index, parnames, pos = ifelse(position == "left", 4, 2)) values.upper.nonInf <- ifelse(values["Upper", ] == Inf, 1.5 * xlim[[2]], values["Upper", ]) # Suppress warnings for non-existing arrow lengths - suppressWarnings(arrows(as.numeric(values["Lower", ]), parname_index, - as.numeric(values.upper.nonInf), parname_index, + suppressWarnings(arrows(as.numeric(values["Lower", ]), parname_index, + as.numeric(values.upper.nonInf), parname_index, code = 3, angle = 90, length = 0.05)) } par(oldpars) diff --git a/R/mkinpredict.R b/R/mkinpredict.R index 5f11c35a..994af703 100644 --- a/R/mkinpredict.R +++ b/R/mkinpredict.R @@ -17,8 +17,8 @@ # You should have received a copy of the GNU General Public License along with # this program. If not, see -mkinpredict <- function(mkinmod, odeparms, odeini, - outtimes, solution_type = "deSolve", +mkinpredict <- function(mkinmod, odeparms, odeini, + outtimes, solution_type = "deSolve", use_compiled = "auto", method.ode = "lsoda", atol = 1e-8, rtol = 1e-10, map_output = TRUE, ...) { @@ -40,12 +40,12 @@ mkinpredict <- function(mkinmod, odeparms, odeini, # Create a function calculating the differentials specified by the model # if necessary if (solution_type == "analytical") { - parent.type = names(mkinmod$map[[1]])[1] + parent.type = names(mkinmod$map[[1]])[1] parent.name = names(mkinmod$diffs)[[1]] o <- switch(parent.type, - SFO = SFO.solution(outtimes, + SFO = SFO.solution(outtimes, evalparse(parent.name), - ifelse(mkinmod$use_of_ff == "min", + ifelse(mkinmod$use_of_ff == "min", evalparse(paste("k", parent.name, "sink", sep="_")), evalparse(paste("k", parent.name, sep="_")))), FOMC = FOMC.solution(outtimes, @@ -53,7 +53,7 @@ mkinpredict <- function(mkinmod, odeparms, odeini, evalparse("alpha"), evalparse("beta")), IORE = IORE.solution(outtimes, evalparse(parent.name), - ifelse(mkinmod$use_of_ff == "min", + ifelse(mkinmod$use_of_ff == "min", evalparse(paste("k__iore", parent.name, "sink", sep="_")), evalparse(paste("k__iore", parent.name, sep="_"))), evalparse("N_parent")), @@ -75,18 +75,18 @@ mkinpredict <- function(mkinmod, odeparms, odeini, names(out) <- c("time", sub("_free", "", parent.name)) } if (solution_type == "eigen") { - coefmat.num <- matrix(sapply(as.vector(mkinmod$coefmat), evalparse), + coefmat.num <- matrix(sapply(as.vector(mkinmod$coefmat), evalparse), nrow = length(mod_vars)) e <- eigen(coefmat.num) c <- solve(e$vectors, odeini) f.out <- function(t) { e$vectors %*% diag(exp(e$values * t), nrow=length(mod_vars)) %*% c } - o <- matrix(mapply(f.out, outtimes), + o <- matrix(mapply(f.out, outtimes), nrow = length(mod_vars), ncol = length(outtimes)) out <- data.frame(outtimes, t(o)) names(out) <- c("time", mod_vars) - } + } if (solution_type == "deSolve") { if (!is.null(mkinmod$cf) & use_compiled[1] != FALSE) { out <- ode( @@ -108,7 +108,7 @@ mkinpredict <- function(mkinmod, odeparms, odeini, diffs <- vector() for (box in names(mkinmod$diffs)) { - diffname <- paste("d", box, sep="_") + diffname <- paste("d", box, sep="_") diffs[diffname] <- with(as.list(c(time, state, parms)), eval(parse(text=mkinmod$diffs[[box]]))) } @@ -117,14 +117,14 @@ mkinpredict <- function(mkinmod, odeparms, odeini, out <- ode( y = odeini, times = outtimes, - func = mkindiff, + func = mkindiff, parms = odeparms, method = method.ode, atol = atol, rtol = rtol, ... ) - } + } if (sum(is.na(out)) > 0) { stop("Differential equations were not integrated for all output times because\n", "NaN values occurred in output from ode()") @@ -140,7 +140,7 @@ mkinpredict <- function(mkinmod, odeparms, odeini, out_mapped[var] <- rowSums(out[, mkinmod$map[[var]]]) } } - return(out_mapped) + return(out_mapped) } else { return(out) } diff --git a/R/mkinresplot.R b/R/mkinresplot.R index 82ffd2cb..3650ef4b 100644 --- a/R/mkinresplot.R +++ b/R/mkinresplot.R @@ -17,16 +17,16 @@ # this program. If not, see if(getRversion() >= '2.15.1') utils::globalVariables(c("variable", "residual")) -mkinresplot <- function (object, +mkinresplot <- function (object, obs_vars = names(object$mkinmod$map), xlim = c(0, 1.1 * max(object$data$time)), xlab = "Time", ylab = "Residual", - maxabs = "auto", legend= TRUE, lpos = "topright", ...) + maxabs = "auto", legend= TRUE, lpos = "topright", ...) { obs_vars_all <- as.character(unique(object$data$variable)) if (length(obs_vars) > 0){ - obs_vars <- intersect(obs_vars_all, obs_vars) + obs_vars <- intersect(obs_vars_all, obs_vars) } else obs_vars <- obs_vars_all residuals <- subset(object$data, variable %in% obs_vars, residual) @@ -37,7 +37,7 @@ mkinresplot <- function (object, names(col_obs) <- names(pch_obs) <- obs_vars plot(0, type = "n", - xlab = xlab, ylab = ylab, + xlab = xlab, ylab = ylab, xlim = xlim, ylim = c(-1.2 * maxabs, 1.2 * maxabs), ...) @@ -49,7 +49,7 @@ mkinresplot <- function (object, abline(h = 0, lty = 2) if (legend == TRUE) { - legend(lpos, inset = c(0.05, 0.05), legend = obs_vars, + legend(lpos, inset = c(0.05, 0.05), legend = obs_vars, col = col_obs[obs_vars], pch = pch_obs[obs_vars]) } } diff --git a/R/mmkin.R b/R/mmkin.R index be4526af..63542f4f 100644 --- a/R/mmkin.R +++ b/R/mmkin.R @@ -19,10 +19,10 @@ # You should have received a copy of the GNU General Public License along with # this program. If not, see -mmkin <- function(models = c("SFO", "FOMC", "DFOP"), datasets, - cores = round(detectCores()/2), cluster = NULL, ...) +mmkin <- function(models = c("SFO", "FOMC", "DFOP"), datasets, + cores = round(detectCores()/2), cluster = NULL, ...) { - parent_models_available = c("SFO", "FOMC", "DFOP", "HS", "SFORB", "IORE") + parent_models_available = c("SFO", "FOMC", "DFOP", "HS", "SFORB", "IORE") n.m <- length(models) n.d <- length(datasets) n.fits <- n.m * n.d @@ -32,10 +32,10 @@ mmkin <- function(models = c("SFO", "FOMC", "DFOP"), datasets, if (!all(sapply(models, function(x) inherits(x, "mkinmod")))) { if (!all(models %in% parent_models_available)) { stop("Please supply models as a list of mkinmod objects or a vector combined of\n ", - paste(parent_models_available, collapse = ", ")) + paste(parent_models_available, collapse = ", ")) } else { names(models) <- models - } + } } else { if (is.null(names(models))) names(models) <- as.character(1:n.m) } diff --git a/R/plot.mkinfit.R b/R/plot.mkinfit.R index 00598737..3df6c81c 100644 --- a/R/plot.mkinfit.R +++ b/R/plot.mkinfit.R @@ -20,15 +20,15 @@ if(getRversion() >= '2.15.1') utils::globalVariables(c("type", "variable", "obse plot.mkinfit <- function(x, fit = x, obs_vars = names(fit$mkinmod$map), xlab = "Time", ylab = "Observed", - xlim = range(fit$data$time), + xlim = range(fit$data$time), ylim = "default", col_obs = 1:length(obs_vars), - pch_obs = col_obs, + pch_obs = col_obs, lty_obs = rep(1, length(obs_vars)), - add = FALSE, legend = !add, + add = FALSE, legend = !add, show_residuals = FALSE, maxabs = "auto", sep_obs = FALSE, rel.height.middle = 0.9, - lpos = "topright", inset = c(0.05, 0.05), + lpos = "topright", inset = c(0.05, 0.05), show_errmin = FALSE, errmin_digits = 3, ...) { if (add && show_residuals) stop("If adding to an existing plot we can not show residuals") @@ -53,12 +53,12 @@ plot.mkinfit <- function(x, fit = x, rownames(subset(fit$fixed, type == "deparm"))) odeparms <- parms.all[odenames] - out <- try(mkinpredict(fit$mkinmod, odeparms, odeini, outtimes, + out <- try(mkinpredict(fit$mkinmod, odeparms, odeini, outtimes, solution_type = solution_type, atol = fit$atol, rtol = fit$rtol), silent = TRUE) if (inherits(out, "try-error")) { - out <- mkinpredict(fit$mkinmod, odeparms, odeini, outtimes, + out <- mkinpredict(fit$mkinmod, odeparms, odeini, outtimes, solution_type = solution_type, atol = fit$atol, rtol = fit$rtol, use_compiled = FALSE) } @@ -85,7 +85,7 @@ plot.mkinfit <- function(x, fit = x, # and the middle plots (if n_plot_rows >2) are smaller by rel.height.middle rel.heights <- if (n_plot_rows > 2) c(1, rep(rel.height.middle, n_plot_rows - 2), 1) else rep(1, n_plot_rows) - layout_matrix = matrix(1:n_plots, + layout_matrix = matrix(1:n_plots, n_plot_rows, n_plot_cols, byrow = TRUE) layout(layout_matrix, heights = rel.heights) } else { # else show residuals in the lower third to keep compatibility @@ -104,8 +104,8 @@ plot.mkinfit <- function(x, fit = x, # Set ylim to sensible default, or to the specified value if (ylim[[1]] == "default") { - ylim_row = c(0, max(c(subset(fit$data, variable %in% row_obs_vars)$observed, - subset(fit$data, variable %in% row_obs_vars)$fitted), + ylim_row = c(0, max(c(subset(fit$data, variable %in% row_obs_vars)$observed, + subset(fit$data, variable %in% row_obs_vars)$fitted), na.rm = TRUE)) } else { ylim_row = ylim @@ -115,12 +115,12 @@ plot.mkinfit <- function(x, fit = x, # Margins for top row of plots when we have more than one row # Reduce bottom margin by 2.1 - hides x axis legend if (plot_row == 1 & n_plot_rows > 1) { - par(mar = c(3.0, 4.1, 4.1, 2.1)) + par(mar = c(3.0, 4.1, 4.1, 2.1)) } # Margins for middle rows of plots, if any if (plot_row > 1 & plot_row < n_plot_rows) { - # Reduce top margin by 2 after the first plot as we have no main title, + # Reduce top margin by 2 after the first plot as we have no main title, # reduced plot height, therefore we need rel.height.middle in the layout par(mar = c(3.0, 4.1, 2.1, 2.1)) } @@ -134,14 +134,14 @@ plot.mkinfit <- function(x, fit = x, # Set up the main plot if not to be added to an existing plot if (add == FALSE) { - plot(0, type="n", + plot(0, type="n", xlim = xlim, ylim = ylim_row, xlab = xlab, ylab = ylab, ...) } # Plot the data for (obs_var in row_obs_vars) { - points(subset(fit$data, variable == obs_var, c(time, observed)), + points(subset(fit$data, variable == obs_var, c(time, observed)), pch = pch_obs[obs_var], col = col_obs[obs_var]) } @@ -181,8 +181,8 @@ plot.mkinfit <- function(x, fit = x, residuals <- subset(fit$data, variable %in% row_obs_vars, residual) if (maxabs == "auto") maxabs = max(abs(residuals), na.rm = TRUE) if (!sep_obs) par(mar = c(5, 4, 0, 2) + 0.1) - plot(0, type="n", - xlim = xlim, + plot(0, type="n", + xlim = xlim, ylim = c(-1.2 * maxabs, 1.2 * maxabs), xlab = xlab, ylab = "Residuals") for(obs_var in row_obs_vars){ @@ -197,6 +197,6 @@ plot.mkinfit <- function(x, fit = x, # Convenience function for switching on some features of mkinfit # that have not been made the default to keep compatibility plot_sep <- function(fit, sep_obs = TRUE, show_residuals = TRUE, show_errmin = TRUE, ...) { - plot.mkinfit(fit, sep_obs = TRUE, show_residuals = TRUE, + plot.mkinfit(fit, sep_obs = TRUE, show_residuals = TRUE, show_errmin = TRUE, ...) } diff --git a/R/plot.mmkin.R b/R/plot.mmkin.R index 7de91e3e..562bbb71 100644 --- a/R/plot.mmkin.R +++ b/R/plot.mmkin.R @@ -16,12 +16,12 @@ # You should have received a copy of the GNU General Public License along with # this program. If not, see -plot.mmkin <- function(x, main = "auto", legends = 1, errmin_var = "All data", errmin_digits = 3, +plot.mmkin <- function(x, main = "auto", legends = 1, errmin_var = "All data", errmin_digits = 3, cex = 0.7, rel.height.middle = 0.9, ...) { n.m <- nrow(x) n.d <- ncol(x) - # We can handle either a row (different models, same dataset) + # We can handle either a row (different models, same dataset) # or a column (same model, different datasets) if (n.m > 1 & n.d > 1) stop("Please select fits either for one model or for one dataset") if (n.m == 1 & n.d == 1) loop_over = "none" @@ -53,12 +53,12 @@ plot.mmkin <- function(x, main = "auto", legends = 1, errmin_var = "All data", e # Margins for top row of plots when we have more than one row # Reduce bottom margin by 2.1 - hides x axis legend if (i.fit == 1 & n.fits > 1) { - par(mar = c(3.0, 4.1, 4.1, 2.1)) + par(mar = c(3.0, 4.1, 4.1, 2.1)) } # Margins for middle rows of plots, if any if (i.fit > 1 & i.fit < n.fits) { - # Reduce top margin by 2 after the first plot as we have no main title, + # Reduce top margin by 2 after the first plot as we have no main title, # reduced plot height, therefore we need rel.height.middle in the layout par(mar = c(3.0, 4.1, 2.1, 2.1)) } @@ -77,7 +77,7 @@ plot.mmkin <- function(x, main = "auto", legends = 1, errmin_var = "All data", e fit_name <- switch(loop_over, models = rownames(x)[i.fit], datasets = colnames(x)[i.fit], - none = "") + none = "") chi2 <- paste0(signif(100 * mkinerrmin(fit)[errmin_var, "err.min"], errmin_digits), "%") mtext(bquote(.(fit_name) ~ chi^2 ~ "error level" == .(chi2)), cex = cex, line = 0.4) diff --git a/R/transform_odeparms.R b/R/transform_odeparms.R index 0946ff14..c871c52a 100644 --- a/R/transform_odeparms.R +++ b/R/transform_odeparms.R @@ -17,10 +17,10 @@ # this program. If not, see transform_odeparms <- function(parms, mkinmod, - transform_rates = TRUE, - transform_fractions = TRUE) + transform_rates = TRUE, + transform_fractions = TRUE) { - # We need the model specification for the names of the model + # We need the model specification for the names of the model # variables and the information on the sink spec = mkinmod$spec @@ -79,7 +79,7 @@ transform_odeparms <- function(parms, mkinmod, } else { transparms[pname] <- parms[pname] } - } + } } # DFOP parameter g is treated as a fraction @@ -99,7 +99,7 @@ backtransform_odeparms <- function(transparms, mkinmod, transform_rates = TRUE, transform_fractions = TRUE) { - # We need the model specification for the names of the model + # We need the model specification for the names of the model # variables and the information on the sink spec = mkinmod$spec -- cgit v1.2.1