# $Id$
# Some of the CAKE R modules are based on mkin, 
# Developed by Hybrid Intelligence (formerly Tessella), part of Capgemini Engineering,
# for Syngenta: Copyright (C) 2011-2022 Syngenta
# Tessella Project Reference: 6245, 7247, 8361, 7414, 10091

#    The CAKE R modules are free software: you can redistribute it and/or modify
#    it under the terms of the GNU General Public License as published by
#    the Free Software Foundation, either version 3 of the License, or
#    (at your option) any later version.
# 
#    This program is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#    GNU General Public License for more details.
# 
#    You should have received a copy of the GNU General Public License
#    along with this program.  If not, see <http://www.gnu.org/licenses/>.

# Shifts parameters slightly away from boundaries specified in "lower" and 
# "upper" (to avoid computational issues after parameter transforms in modFit).
ShiftAwayFromBoundaries <- function(parameters, lower, upper) {
    parametersOnLowerBound = which(parameters == lower)
    parameters[parametersOnLowerBound] <- parameters[parametersOnLowerBound] * (1 + .Machine$double.eps) + .Machine$double.xmin

    parametersOnUpperBound = which(parameters == upper)
    parameters[parametersOnUpperBound] <- parameters[parametersOnUpperBound] * (1 - .Machine$double.neg.eps) - .Machine$double.xmin

    return(parameters)
}

# Adjusts stated initial values to put into the ODE solver.
#
# odeini: The initial values to adjust (in the form that would be fed into the ode function).
# cake.model: The expression of the model that we are solving.
# odeparms: The parameters for the ODE (in the form that would be fed into the ode function).
#
# Returns: Adjusted initial values.
AdjustOdeInitialValues <- function(odeini, cake.model, odeparms) {
    odeini.names <- names(odeini)

    for (ini.name in odeini.names) {
        # For DFOP metabolites in two compartments, need to calculate some initial conditions for the ODEs.
        if (!(ini.name %in% names(cake.model$diffs))) {
            subcompartment1.name <- paste(ini.name, "1", sep = "_")
            subcompartment2.name <- paste(ini.name, "2", sep = "_")

            if (subcompartment1.name %in% names(cake.model$diffs) && subcompartment2.name %in% names(cake.model$diffs)) {
                g.parameter.name = paste("g", ini.name, sep = "_")

                odeini[[subcompartment1.name]] <- odeini[[ini.name]] * odeparms[[g.parameter.name]]
                odeini[[subcompartment2.name]] <- odeini[[ini.name]] * (1 - odeparms[[g.parameter.name]])
            }
        }
    }

    # It is important that these parameters are stated in the same order as the differential equations.
    return(odeini[names(cake.model$diffs)])
}

# Post-processes the output from the ODE solver (or analytical process), including recombination of sub-compartments.
#
# odeoutput: The output of the ODE solver.
# cake.model: The expression of the model that we are solving.
# atol: The tolerance to which the solution has been calculated.
#
# Returns: Post-processed/transformed ODE output.
PostProcessOdeOutput <- function(odeoutput, cake.model, atol) {
    out_transformed <- data.frame(time = odeoutput[, "time"])

    # Replace values that are incalculably small with 0.
    for (col.name in colnames(odeoutput)) {
        if (col.name == "time") {
            next
        }

        # If we have non-NaN, positive outputs...
        if (length(odeoutput[, col.name][!is.nan(odeoutput[, col.name]) & odeoutput[, col.name] > 0]) > 0) {
            # ...then replace the NaN outputs.
            odeoutput[, col.name][is.nan(odeoutput[, col.name])] <- 0
        }

        # Round outputs smaller than the used tolerance down to 0.
        odeoutput[, col.name][odeoutput[, col.name] < atol] <- 0
    }

    # Re-combine sub-compartments (if required)
    for (compartment.name in names(cake.model$map)) {
        if (length(cake.model$map[[compartment.name]]) == 1) {
            out_transformed[compartment.name] <- odeoutput[, compartment.name]
        } else {
            out_transformed[compartment.name] <- rowSums(odeoutput[, cake.model$map[[compartment.name]]])
        }
    }

    return(out_transformed)
}

# Reorganises data in a wide format to a long format.
#
# wide_data: The data in wide format.
# time: The name of the time variable in wide_data (default "t").
#
# Returns: Reorganised data.
wide_to_long <- function(wide_data, time = "t") {
    colnames <- names(wide_data)

    if (!(time %in% colnames)) {
        stop("The data in wide format have to contain a variable named ", time, ".")
    }

    vars <- subset(colnames, colnames != time)
    n <- length(colnames) - 1
    long_data <- data.frame(name = rep(vars, each = length(wide_data[[time]])),
                          time = as.numeric(rep(wide_data[[time]], n)), value = as.numeric(unlist(wide_data[vars])),
                          row.names = NULL)

    return(long_data)
}

RunFitStep <- function(cost, costForExtraSolver, useExtraSolver, parameters, lower, upper, control) {
    if (useExtraSolver) {
        a <- try(fit <- solnp(parameters, fun = costForExtraSolver, LB = lower, UB = upper, control = control), silent = TRUE)

        fitted_with_extra_solver <- TRUE
        if (class(a) == "try-error") {
            cat('Extra solver failed, trying PORT')
            ## now using submethod already
            a <- try(fit <- modFit(cost, parameters, lower = lower, upper = upper, method = 'Port', control = control))
            fitted_with_extra_solver <- FALSE
            if (class(a) == "try-error") {
                cat('PORT failed, trying L-BFGS-B')
                fit <- modFit(cost, parameters, lower = lower, upper = upper, method = 'L-BFGS-B', control = control)
            }
        }
    } else {
        # modFit parameter transformations can explode if you put in parameters that are equal to a bound, so we move them away by a tiny amount.
        all.optim <- ShiftAwayFromBoundaries(parameters, lower, upper)

        fit <- modFit(cost, all.optim, lower = lower,
                          upper = upper, control = control)
        fitted_with_extra_solver <- FALSE
    }
    return(list(fit = fit, fitted_with_extra_solver = fitted_with_extra_solver))
}

GetFitValuesAfterExtraSolver <- function(fit, cake_cost) {
    fit$ssr <- fit$values[length(fit$values)]
    fit$residuals <- cake_cost$residual$res
    ## mean square per varaible
    if (class(cake_cost) == "modCost") {
        names(fit$residuals) <- cake_cost$residuals$name
        fit$var_ms <- cake_cost$var$SSR / cake_cost$var$N
        fit$var_ms_unscaled <- cake_cost$var$SSR.unscaled / cake_cost$var$N
        fit$var_ms_unweighted <- cake_cost$var$SSR.unweighted / cake_cost$var$N
        names(fit$var_ms_unweighted) <- names(fit$var_ms_unscaled) <-
                    names(fit$var_ms) <- FF$var$name
    } else fit$var_ms <- fit$var_ms_unweighted <- fit$var_ms_unscaled <- NA

    return(fit)
}

GetOptimiserSpecificSetup <- function(optimiser) {
    switch(optimiser,
        OLS = GetOlsSpecificSetup(),
        IRLS = GetIrlsSpecificSetup(),
        MCMC = GetMcmcSpecificSetup())    
}

GetOptimisationRoutine <- function(optimiser) {
    switch(optimiser,
        OLS = GetOlsOptimisationRoutine(),
        IRLS = GetIrlsOptimisationRoutine(),
        MCMC = GetMcmcOptimisationRoutine())
}

GetOptimiserSpecificWrapUp <- function(optimiser) {
    switch(optimiser,
        OLS = GetOlsSpecificWrapUp(),
        IRLS = GetIrlsSpecificWrapUp(),
        MCMC = GetMcmcSpecificWrapUp())
}