path: root/CakeOlsFit.R

# Originally: mkinfit.R 87 2010-12-09 07:31:59Z jranke $

# Based on code in mkinfit
# Portions Johannes Ranke 2010
# Contact: mkin-devel@lists.berlios.de
# The summary function is an adapted and extended version of summary.modFit
# from the FME package, v 1.1 by Soetart and Petzoldt, which was in turn
# inspired by summary.nls.lm

# Modifications 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/>.


# Performs an ordinary least squares fit on a given CAKE model.
#
# 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).
# atol: The tolerance to apply to the ODE solver.
# dfopDtMaxIter: The maximum number of iterations to apply to DFOP DT calculation.
CakeOlsFit <- function(cake.model,
                       observed,
                       parms.ini,
                       state.ini,
                       lower = 0,
                       upper = Inf,
                       fixed_parms = NULL,
                       fixed_initials = names(cake.model$diffs)[-1],
                       quiet = FALSE,
                       atol = 1e-6,
                       dfopDtMaxIter = 10000,
                       control = list(),
                       useExtraSolver = FALSE) {
    fit <- CakeFit("OLS",
                    cake.model,
                    observed,
                    parms.ini,
                    state.ini,
                    lower,
                    upper,
                    fixed_parms,
                    fixed_initials,
                    quiet,
                    atol = atol,
                    dfopDtMaxIter = dfopDtMaxIter,
                    control = control,
                    useExtraSolver = useExtraSolver)

    return(fit)
}

GetOlsSpecificSetup <- function() {
    return(function(
           cake.model,
           state.ini.optim,
           state.ini.optim.boxnames,
           state.ini.fixed,
           parms.fixed,
           observed,
           mkindiff,
           quiet,
           atol,
           solution,
           ...) {
        costFunctions <- CakeInternalCostFunctions(cake.model, state.ini.optim, state.ini.optim.boxnames, state.ini.fixed,
                                             parms.fixed, observed, mkindiff = mkindiff, quiet, atol = atol, solution = solution, err = NULL)
        return(list(costFunctions = costFunctions, costWithStatus = NULL))
    })
}

GetOlsOptimisationRoutine <- function() {
    return(function(costFunctions, costForExtraSolver, useExtraSolver, parms, lower, upper, control, ...) {
        fitStepResult <- RunFitStep(costFunctions$cost, costForExtraSolver, useExtraSolver, parms, lower, upper, control)
        fit <- fitStepResult$fit
        fitted_with_extra_solver <- fitStepResult$fitted_with_extra_solver

        if (fitted_with_extra_solver) {
            fit <- GetFitValuesAfterExtraSolver(fit, FF)

            np <- length(parms)
            fit$rank <- np
            fit$df.residual <- length(fit$residuals) - fit$rank

            # solnp can return an incorrect Hessian, so we use another fitting method at the optimised point to determine the Hessian
            fit <- modFit(costFunctions$cost, fit$par, lower = lower, upper = upper, method = 'L-BFGS-B', control = list())
        }
        return(list(fit = fit, fitted_with_extra_solver = fitted_with_extra_solver, res = NULL))
    })
}

GetOlsSpecificWrapUp <- function() {
    return(function(fit, ...) {
        args <- list(...)
        parms.fixed <- args$parms.fixed
        observed <- args$observed

        parms.all <- c(fit$par, parms.fixed)

        data <- observed
        data$err <- rep(NA, length(data$time))

        class(fit) <- c("CakeFit", "mkinfit", "modFit")

        return(list(fit = fit, parms.all = parms.all, data = data))
    })
}