diff options
Diffstat (limited to 'CakeFit.R')
| -rwxr-xr-x | CakeFit.R | 205 |
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 |