summaryrefslogtreecommitdiff
path: root/CakeIrlsFit.R
diff options
context:
space:
mode:
Diffstat (limited to 'CakeIrlsFit.R')
-rwxr-xr-x[-rw-r--r--]CakeIrlsFit.R394
1 files changed, 140 insertions, 254 deletions
diff --git a/CakeIrlsFit.R b/CakeIrlsFit.R
index f1629f2..6a50621 100644..100755
--- a/CakeIrlsFit.R
+++ b/CakeIrlsFit.R
@@ -1,7 +1,8 @@
#$Id$
#
# Some of the CAKE R modules are based on mkin
-# 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
@@ -34,269 +35,154 @@
# dfopDtMaxIter: The maximum number of iterations to apply to DFOP DT calculation.
# control: ...
# useExtraSolver: Whether to use the extra solver for this fit.
-CakeIrlsFit <- function (cake.model,
- observed,
- parms.ini = rep(0.1, length(cake.model$parms)),
- state.ini = c(100, rep(0, length(cake.model$diffs) - 1)),
- 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,
- ...)
-{
- NAind <-which(is.na(observed$value))
- mod_vars <- names(cake.model$diffs)
- observed <- subset(observed, name %in% names(cake.model$map))
- ERR <- rep(1,nrow(observed))
- observed <- cbind(observed,err=ERR)
- fitted_with_extra_solver <- 0
-
- obs_vars <- unique(as.character(observed$name))
-
- if (is.null(names(parms.ini))) {
- names(parms.ini) <- cake.model$parms
- }
-
- 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)))
- }
-
- if (is.null(names(state.ini))) {
- names(state.ini) <- mod_vars
- }
-
- 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 = "_")
- }
-
- costFunctions <- CakeInternalCostFunctions(cake.model, state.ini.optim, state.ini.optim.boxnames,
- state.ini.fixed, parms.fixed, observed, mkindiff, quiet, atol=atol)
-
- ############### Iteratively Reweighted Least Squares#############
- ## Start with no weighting
-
- ## a prefitting step since this is usually the most effective method
- if(useExtraSolver)
- {
- pnames=names(c(state.ini.optim, parms.optim))
- fn <- function(P){
- names(P) <- pnames
- FF<<-costFunctions$cost(P)
- return(FF$model)}
- a <- try(fit <- solnp(c(state.ini.optim, parms.optim),fun=fn,LB=lower,UB=upper,control=control),silent=TRUE)
+CakeIrlsFit <- 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("IRLS",
+ cake.model,
+ observed,
+ parms.ini,
+ state.ini,
+ lower,
+ upper,
+ fixed_parms,
+ fixed_initials,
+ quiet,
+ atol = atol,
+ dfopDtMaxIter = dfopDtMaxIter,
+ control = control,
+ useExtraSolver = useExtraSolver)
+
+ return(fit)
+}
+
+GetIrlsSpecificSetup <- function() {
+ return(function(
+ cake.model,
+ state.ini.optim,
+ state.ini.optim.boxnames,
+ state.ini.fixed,
+ parms.fixed,
+ observed,
+ mkindiff,
+ quiet,
+ atol,
+ solution,
+ ...) {
+ control <- list(...)$control
+ # Get the CAKE cost functions
+ costFunctions <- CakeInternalCostFunctions(cake.model, state.ini.optim, state.ini.optim.boxnames, state.ini.fixed,
+ parms.fixed, observed, mkindiff = mkindiff, quiet, atol = atol, solution = solution)
- fitted_with_extra_solver <- 1
- if(class(a) == "try-error")
- {
- print('solnp fails, try PORT or other algorithm by users choice, might take longer time. Do something else!')
- warning('solnp fails, switch to PORT or other algorithm by users choice')
- ## now using submethod already
- a <- try(fit <- modFit(costFunctions$cost, c(state.ini.optim, parms.optim), lower = lower, upper = upper, method='Port',control=control))
- fitted_with_extra_solver <- 0
- if(class(a) == "try-error")
- {
- fit <- modFit(costFunctions$cost, c(state.ini.optim, parms.optim), lower = lower, upper = upper, method='L-BFGS-B',control=control)
+ if (length(control) == 0) {
+ irls.control <- list(maxIter = 10, tol = 1e-05)
+ control <- list(irls.control = irls.control)
+ } else {
+ if (is.null(control$irls.control)) {
+ irls.control <- list(maxIter = 10, tol = 1e-05)
+ control <- list(irls.control = irls.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(c(state.ini.optim, parms.optim), lower, upper)
-
- fit <- modFit(costFunctions$cost, all.optim, lower = lower,
- upper = upper,control=control,...)
- }
-
- if(length(control)==0)
- {
- irls.control <- list(maxIter=10,tol=1e-05)
- control <- list(irls.control=irls.control)
- }else{
- if(is.null(control$irls.control))
- {
- irls.control <- list(maxIter=10,tol=1e-05)
- control <- list(irls.control=irls.control)
- }
- }
-
- irls.control <- control$irls.control
- maxIter <- irls.control$maxIter
- tol <- irls.control$tol
- ####
-
- if(length(cake.model$map)==1){
- ## there is only one parent just do one iteration:
- maxIter <- 0
-
- if(fitted_with_extra_solver==1)## managed to fit with extra solver
- {
- fit$ssr <- fit$values[length(fit$values)]
- fit$residuals <-FF$residual$res
- ## mean square per variable
- if (class(FF) == "modCost") {
- names(fit$residuals) <- FF$residuals$name
- fit$var_ms <- FF$var$SSR/FF$var$N
- fit$var_ms_unscaled <- FF$var$SSR.unscaled/FF$var$N
- fit$var_ms_unweighted <- FF$var$SSR.unweighted/FF$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
+ irls.control <- control$irls.control
+ maxIter <- irls.control$maxIter
+ tol <- irls.control$tol
+
+ return(list(costFunctions = costFunctions, control = control, maxIter = maxIter, tol = tol, costWithStatus = NULL))
+ })
+}
+
+GetIrlsOptimisationRoutine <- function() {
+ return(function(costFunctions, costForExtraSolver, useExtraSolver, parms, lower, upper, control, ...) {
+ irlsArgs <- list(...)
+ cake.model <- irlsArgs$cake.model
+ tol <- irlsArgs$tol
+ fitted_with_extra_solver <- irlsArgs$fitted_with_extra_solver
+ observed <- irlsArgs$observed
+ maxIter <- irlsArgs$maxIter
+
+ pnames = names(parms)
+
+ fitStepResult <- RunFitStep(costFunctions$cost, costForExtraSolver, useExtraSolver, parms, lower, upper, control)
+ fit <- fitStepResult$fit
+ fitted_with_extra_solver <- fitStepResult$fitted_with_extra_solver
+
+ if (length(cake.model$map) == 1) {
+ ## there is only one parent then don't do any re-weighting
+ maxIter <- 0
}
- err1 <- sqrt(fit$var_ms_unweighted)
- ERR <- err1[as.character(observed$name)]
- observed$err <-ERR
- }
-
- niter <- 1
- ## insure one IRLS iteration
- diffsigma <- 100
- olderr <- rep(1,length(cake.model$map))
- while(diffsigma>tol & niter<=maxIter)
- {
- if(fitted_with_extra_solver==1 && useExtraSolver)## managed to fit with extra solver
- {
- fit$ssr <- fit$values[length(fit$values)]
- fit$residuals <-FF$residual$res
- ## mean square per variable
- if (class(FF) == "modCost") {
- names(fit$residuals) <- FF$residuals$name
- fit$var_ms <- FF$var$SSR/FF$var$N
- fit$var_ms_unscaled <- FF$var$SSR.unscaled/FF$var$N
- fit$var_ms_unweighted <- FF$var$SSR.unweighted/FF$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
+ niter <- 1
+ ## ensure one IRLS iteration
+ diffsigma <- 100
+ olderr <- rep(1, length(cake.model$map))
+ while (diffsigma > tol & niter <= maxIter) {
+ # Read info from FF into fit if extra solver was used and set observed$err
+ if (fitted_with_extra_solver) {
+ fit <- GetFitValuesAfterExtraSolver(fit, FF)
+ }
+ err <- sqrt(fit$var_ms_unweighted)
+ ERR <- err[as.character(observed$name)]
+ costFunctions$set.error(ERR)
+
+ diffsigma <- sum((err - olderr) ^ 2)
+ cat("IRLS iteration at", niter, "; Diff in error variance ", diffsigma, "\n")
+ olderr <- err
+
+ fitStepResult <- RunFitStep(costFunctions$cost, costForExtraSolver, useExtraSolver, fit$par, lower, upper, control)
+
+ fit <- fitStepResult$fit
+ fitted_with_extra_solver <- fitStepResult$fitted_with_extra_solver
+
+ niter <- niter + 1
+
+ # If not converged, reweight and fit
}
- err <- sqrt(fit$var_ms_unweighted)
- ERR <- err[as.character(observed$name)]
- costFunctions$set.error(ERR)
- diffsigma <- sum((err-olderr)^2)
- cat("IRLS iteration at",niter, "; Diff in error variance ", diffsigma,"\n")
- olderr <- err
-
- if(useExtraSolver)
- {
- fitted_with_extra_solver <- 1
- a <- try(fit <- solnp(fit$par,fun=fn,LB=lower,UB=upper,control=control),silent=TRUE)
+ # Read info from FF into fit if extra solver was used, set fit$residuals and get correct Hessian as solnp doesn't
+ if (fitted_with_extra_solver) {
+ fit <- GetFitValuesAfterExtraSolver(fit, FF)
- if(class(a) == "try-error")
- {
- fitted_with_extra_solver <- 0
- print('solnp fails during IRLS iteration, try PORT or other algorithm by users choice.This may takes a while. Do something else!') ## NOTE: because in kingui we switch off the warnings, we need to print out the message instead.
- warning('solnp fails during IRLS iteration, switch to PORT or other algorithm by users choice')
-
- fit <- modFit(costFunctions$cost, fit$par, lower = lower, upper = upper, method='Port',control=list())
- }
- }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.
- fit$par <- ShiftAwayFromBoundaries(fit$par, lower, upper)
-
- fit <- modFit(costFunctions$cost, fit$par, lower = lower, upper = upper, control=control, ...)
+ 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())
}
-
- niter <- niter+1
-
- ### If not converged, reweight and fit
- }
-
- if(fitted_with_extra_solver==1 && useExtraSolver){
- ## solnp used
- optimmethod <- 'solnp'
- fit$ssr <- fit$values[length(fit$values)]
- fit$residuals <-FF$residual$res
- ## mean square per varaible
- if (class(FF) == "modCost") {
- names(fit$residuals) <- FF$residuals$name
- fit$var_ms <- FF$var$SSR/FF$var$N
- fit$var_ms_unscaled <- FF$var$SSR.unscaled/FF$var$N
- fit$var_ms_unweighted <- FF$var$SSR.unweighted/FF$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
- np <- length(c(state.ini.optim, parms.optim))
- 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())
- }
-
- ###########################################
- 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 = "deSolve", atol)
- out_transformed <- PostProcessOdeOutput(out_predicted, cake.model, atol)
-
- predicted_long <- wide_to_long(out_transformed, time = "time")
- fit$predicted <- out_transformed
-
- fit$errmin <- CakeChi2(cake.model, observed, predicted_long, obs_vars, parms.optim, state.ini.optim, state.ini, parms.ini, fit$fixed)
- parms.all <- c(fit$par, parms.fixed, state.ini)
- fit$penalties <- CakePenaltiesLong(parms.all, out_transformed, observed)
- 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)
- }
+ err1 <- sqrt(fit$var_ms_unweighted)
+ ERR <- err1[as.character(observed$name)]
+ observed$err <- ERR
- fit$ioreRepDT <- CakeIORERepresentativeDT("Parent", parms.all)
- fit$fomcRepDT <- CakeFOMCBackCalculatedDT(parms.all)
-
- data <- merge(observed, 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
-
- class(fit) <- c("CakeFit", "mkinfit", "modFit")
- return(fit)
+ return(list(fit = fit, fitted_with_extra_solver = fitStepResult$fitted_with_extra_solver, observed = observed, res = NULL))
+ })
}
+
+GetIrlsSpecificWrapUp <- function() {
+ return(function(fit, ...) {
+ args <- list(...)
+ parms.fixed <- args$parms.fixed
+ state.ini <- args$state.ini
+ observed <- args$observed
+
+ parms.all <- c(fit$par, parms.fixed, state.ini)
+
+ data <- observed
+
+ class(fit) <- c("CakeFit", "mkinfit", "modFit")
+
+ return(list(fit = fit, parms.all = parms.all, data = data))
+ })
+} \ No newline at end of file

Contact - Imprint