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