Version 3.5, from a fresh installationv3.5

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