summaryrefslogtreecommitdiff
path: root/CakeOlsFit.R
diff options
context:
space:
mode:
Diffstat (limited to 'CakeOlsFit.R')
-rwxr-xr-x[-rw-r--r--]CakeOlsFit.R185
1 files changed, 69 insertions, 116 deletions
diff --git a/CakeOlsFit.R b/CakeOlsFit.R
index e9b8f3c..66bf9b4 100644..100755
--- a/CakeOlsFit.R
+++ b/CakeOlsFit.R
@@ -7,7 +7,8 @@
# from the FME package, v 1.1 by Soetart and Petzoldt, which was in turn
# inspired by summary.nls.lm
-# Modifications developed by Tessella for Syngenta: Copyright (C) 2011-2020 Syngenta
+# 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
@@ -39,8 +40,8 @@
# dfopDtMaxIter: The maximum number of iterations to apply to DFOP DT calculation.
CakeOlsFit <- function(cake.model,
observed,
- parms.ini = rep(0.1, length(cake.model$parms)),
- state.ini = c(100, rep(0, length(cake.model$diffs) - 1)),
+ parms.ini,
+ state.ini,
lower = 0,
upper = Inf,
fixed_parms = NULL,
@@ -48,126 +49,78 @@ CakeOlsFit <- function(cake.model,
quiet = FALSE,
atol = 1e-6,
dfopDtMaxIter = 10000,
- ...)
-{
- mod_vars <- names(cake.model$diffs)
- # Subset dataframe with mapped (modelled) variables
- observed <- subset(observed, name %in% names(cake.model$map))
- # Get names of observed variables
- obs_vars <- unique(as.character(observed$name))
-
- # Name the parameters if they are not named yet
- if(is.null(names(parms.ini))) names(parms.ini) <- cake.model$parms
-
- # Name the inital parameter values if they are not named yet
- if(is.null(names(state.ini))) names(state.ini) <- mod_vars
-
- # 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, the spectral decomposition of the matrix (fundamental system)
- # 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) {
- 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)))
- }
- }
-
- 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)
-
- # 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(c(state.ini.optim, parms.optim), lower, upper)
-
- fit <-modFit(costFunctions$cost, all.optim, lower = lower, upper = upper, ...)
-
- # We need to return some more data for summary and plotting
- fit$solution <- solution
-
- if (solution == "deSolve") {
- fit$mkindiff <- mkindiff
- }
-
- # We also need various other information for summary and plotting
- fit$map <- cake.model$map
- fit$diffs <- cake.model$diffs
-
- # Collect initial parameter values in two dataframes
- 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)))
+ 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)
+}
- # Ensure initial state is at time 0
- obstimes <- unique(c(0,observed$time))
+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))
+ })
+}
- out_predicted <- CakeOdeSolve(fit, obstimes, solution = solution, atol)
- out_transformed <- PostProcessOdeOutput(out_predicted, cake.model, atol)
+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
- predicted_long <- wide_to_long(out_transformed, time = "time")
- fit$predicted <- out_transformed
+ if (fitted_with_extra_solver) {
+ fit <- GetFitValuesAfterExtraSolver(fit, FF)
- # 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)
+ np <- length(parms)
+ fit$rank <- np
+ fit$df.residual <- length(fit$residuals) - fit$rank
- # Calculate dissipation times DT50 and DT90 and formation fractions
- parms.all <- c(fit$par, parms.fixed)
- 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))
+ # 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))
+ })
+}
- for (compartment.name in names(cake.model$map)) {
- type <- names(cake.model$map[[compartment.name]])[1]
+GetOlsSpecificWrapUp <- function() {
+ return(function(fit, ...) {
+ args <- list(...)
+ parms.fixed <- args$parms.fixed
+ observed <- args$observed
- fit$distimes[compartment.name, ] <- CakeDT(type,compartment.name,parms.all,dfopDtMaxIter)
- fit$extraDT50[compartment.name, ] <- CakeExtraDT(type, compartment.name, parms.all)
- }
+ parms.all <- c(fit$par, parms.fixed)
- fit$ioreRepDT <- CakeIORERepresentativeDT("Parent", parms.all)
- fit$fomcRepDT <- CakeFOMCBackCalculatedDT(parms.all)
- fit$penalties <- CakePenaltiesLong(parms.all, out_transformed, observed)
+ data <- observed
+ data$err <- rep(NA, length(data$time))
- # Collect observed, predicted and residuals
- data<-observed
- data$err<-rep(NA,length(data$time))
- 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$atol <- atol
- fit$topology <- cake.model$topology
+ class(fit) <- c("CakeFit", "mkinfit", "modFit")
- class(fit) <- c("CakeFit", "mkinfit", "modFit")
- return(fit)
-}
+ return(list(fit = fit, parms.all = parms.all, data = data))
+ })
+} \ No newline at end of file

Contact - Imprint