summaryrefslogtreecommitdiff
path: root/CakeFit.R
diff options
context:
space:
mode:
Diffstat (limited to 'CakeFit.R')
-rwxr-xr-xCakeFit.R205
1 files changed, 205 insertions, 0 deletions
diff --git a/CakeFit.R b/CakeFit.R
new file mode 100755
index 0000000..7e6c15a
--- /dev/null
+++ b/CakeFit.R
@@ -0,0 +1,205 @@
+# Performs the optimisation routine corresponding to the given optimiser selected.
+# Remark: this function was originally based on the "mkinfit" function, version 0.1.
+#
+# optimiser: A character type indicating which optimiser should be used, one of 'OLS', 'IRLS', 'MCMC'
+# cake.model: The model to perform the fit on (as generated by CakeModel.R).
+# observed: Observation data to fit to.
+# parms.ini: Initial values for the parameters being fitted.
+# state.ini: Initial state (i.e. initial values for concentration, the dependent variable being modelled).
+# lower: Lower bounds to apply to parameters.
+# upper: Upper bound to apply to parameters.
+# fixed_parms: A vector of names of parameters that are fixed to their initial values.
+# fixed_initials: A vector of compartments with fixed initial concentrations.
+# quiet: Whether the internal cost functions should execute more quietly than normal (less output).
+# niter: The number of MCMC iterations to apply.
+# atol: The tolerance to apply to the ODE solver.
+# dfopDtMaxIter: The maximum number of iterations to apply to DFOP DT calculation.
+# control: ...
+# useExtraSolver: Whether to use the extra solver for this fit.
+CakeFit <- function(optimiser,
+ cake.model,
+ observed,
+ parms.ini,
+ state.ini,
+ lower = 0,
+ upper = Inf,
+ fixed_parms = NULL,
+ fixed_initials = names(cake.model$diffs)[-1],
+ quiet = FALSE,
+ niter = 1000,
+ verbose = TRUE,
+ seed = NULL,
+ atol = 1e-6,
+ dfopDtMaxIter = 10000,
+ control = list(),
+ useExtraSolver = FALSE) {
+ env = environment()
+
+ # Get model variable names
+ mod_vars <- names(cake.model$diffs)
+
+ # Subset dataframe with mapped (modelled) variables
+ observed <- subset(observed, name %in% names(cake.model$map))
+
+ ERR <- rep(1, nrow(observed))
+ observed <- cbind(observed, err = ERR)
+
+ # Get names of observed variables
+ obs_vars <- unique(as.character(observed$name))
+
+ # Parameters to be optimised
+ parms.fixed <- parms.ini[fixed_parms]
+ optim_parms <- setdiff(names(parms.ini), fixed_parms)
+ parms.optim <- parms.ini[optim_parms]
+
+ state.ini.fixed <- state.ini[fixed_initials]
+ optim_initials <- setdiff(names(state.ini), fixed_initials)
+ state.ini.optim <- state.ini[optim_initials]
+ state.ini.optim.boxnames <- names(state.ini.optim)
+
+ if (length(state.ini.optim) > 0) {
+ names(state.ini.optim) <- paste(names(state.ini.optim), "0", sep = "_")
+ }
+
+ # Decide if the solution of the model can be based on a simple analytical formula
+ # or a numeric ode solver from the deSolve package
+ if (length(cake.model$map) == 1) {
+ solution <- "analytical"
+ } else {
+ solution <- "deSolve"
+ }
+
+ # Create a function calculating the differentials specified by the model
+ # if necessary
+ if (solution == "deSolve") {
+ mkindiff <- function(t, state, parms) {
+ # numerical method, replace with analytical for parent-only fits
+ time <- t
+ diffs <- vector()
+ for (box in mod_vars) {
+ diffname <- paste("d", box, sep = "_")
+ diffs[diffname] <- with(as.list(c(time, state, parms)),
+ eval(parse(text = cake.model$diffs[[box]])))
+ }
+ return(list(c(diffs)))
+ }
+ }
+
+ fitted_with_extra_solver <- FALSE
+
+ optimiserSpecificSetup <- GetOptimiserSpecificSetup(optimiser)
+ optimiserSpecificSetupOutput <- optimiserSpecificSetup(cake.model,
+ state.ini.optim,
+ state.ini.optim.boxnames,
+ state.ini.fixed,
+ parms.fixed,
+ observed,
+ mkindiff = mkindiff,
+ quiet,
+ atol = atol,
+ solution = solution,
+ control = control,
+ seed = seed)
+ list2env(optimiserSpecificSetupOutput, env)
+
+ pnames <- names(c(state.ini.optim, parms.optim))
+ costForExtraSolver <- function(P) {
+ names(P) <- pnames
+ FF <<- costFunctions$cost(P)
+ return(FF$model)
+ }
+
+ optimiserFitRoutine <- GetOptimisationRoutine(optimiser)
+ optimiserFitResult <- optimiserFitRoutine(costFunctions,
+ costForExtraSolver,
+ useExtraSolver,
+ c(state.ini.optim, parms.optim),
+ lower,
+ upper,
+ control,
+ cake.model = cake.model,
+ tol = tol,
+ fitted_with_extra_solver = fitted_with_extra_solver,
+ observed = observed,
+ maxIter = maxIter,
+ costWithStatus = costWithStatus,
+ niter = niter,
+ verbose = verbose)
+ list2env(optimiserFitResult, env)
+
+ optimiserSpecificWrapUp <- GetOptimiserSpecificWrapUp(optimiser)
+ optimiserSpecificWrapUpOutput <- optimiserSpecificWrapUp(fit,
+ parms.fixed = parms.fixed,
+ state.ini = state.ini,
+ observed = observed,
+ res = res,
+ seed = seed,
+ costWithStatus = costWithStatus)
+ list2env(optimiserSpecificWrapUpOutput, env)
+
+ # We need to return some more data for summary and plotting
+ fit$solution <- solution
+
+ if (solution == "deSolve") {
+ fit$mkindiff <- mkindiff
+ }
+
+ fit$map <- cake.model$map
+ fit$diffs <- cake.model$diffs
+ fit$topology <- cake.model$topology
+ fit$start <- data.frame(initial = c(state.ini.optim, parms.optim))
+ fit$start$type <- c(rep("state", length(state.ini.optim)), rep("deparm", length(parms.optim)))
+ fit$start$lower <- lower
+ fit$start$upper <- upper
+
+ fit$fixed <- data.frame(value = c(state.ini.fixed, parms.fixed))
+ fit$fixed$type <- c(rep("state", length(state.ini.fixed)), rep("deparm", length(parms.fixed)))
+
+ # Ensure initial state is at time 0
+ obstimes <- unique(c(0, observed$time))
+
+ out_predicted <- CakeOdeSolve(fit, obstimes, solution = solution, atol)
+ out_transformed <- PostProcessOdeOutput(out_predicted, cake.model, atol)
+
+ predicted_long <- wide_to_long(out_transformed, time = "time")
+ fit$predicted <- out_transformed
+
+ # Calculate chi2 error levels according to FOCUS (2006)
+ fit$errmin <- CakeChi2(cake.model, observed, predicted_long, obs_vars, parms.optim, state.ini.optim, state.ini, parms.ini, fit$fixed)
+
+ # Calculate dissipation times DT50 and DT90 and formation fractions
+ fit$distimes <- data.frame(DT50 = rep(NA, length(obs_vars)), DT90 = rep(NA, length(obs_vars)), row.names = obs_vars)
+ fit$extraDT50 <- data.frame(k1 = rep(NA, length(names(cake.model$map))), k2 = rep(NA, length(names(cake.model$map))), row.names = names(cake.model$map))
+
+ for (compartment.name in names(cake.model$map)) {
+ type <- names(cake.model$map[[compartment.name]])[1]
+
+ fit$distimes[compartment.name,] <- CakeDT(type, compartment.name, parms.all, dfopDtMaxIter)
+ fit$extraDT50[compartment.name,] <- CakeExtraDT(type, compartment.name, parms.all)
+ }
+
+ fit$ioreRepDT <- CakeIORERepresentativeDT("Parent", parms.all)
+ fit$fomcRepDT <- CakeFOMCBackCalculatedDT(parms.all)
+ fit$penalties <- CakePenaltiesLong(parms.all, out_transformed, observed)
+
+ # Collect observed, predicted and residuals
+ data <- merge(data, predicted_long, by = c("time", "name"))
+
+ names(data) <- c("time", "variable", "observed", "err-var", "predicted")
+ data$residual <- data$observed - data$predicted
+ data$variable <- ordered(data$variable, levels = obs_vars)
+ fit$data <- data[order(data$variable, data$time),]
+
+ fit$costFunctions <- costFunctions
+ fit$upper <- upper
+ fit$lower <- lower
+
+ fit$atol <- atol
+
+ if (optimiser == "MCMC") {
+ sq <- data$residual * data$residual
+ fit$ssr <- sum(sq)
+ }
+
+ return(fit)
+} \ No newline at end of file

Contact - Imprint