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