From fd9ea462ba438c023e33902449429eae6d2dc28f Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Wed, 9 Nov 2022 10:08:37 +0100 Subject: Version 3.5, from a fresh installation --- CakeIrlsFit.R | 394 +++++++++++++++++++++------------------------------------- 1 file changed, 140 insertions(+), 254 deletions(-) mode change 100644 => 100755 CakeIrlsFit.R (limited to 'CakeIrlsFit.R') diff --git a/CakeIrlsFit.R b/CakeIrlsFit.R old mode 100644 new mode 100755 index f1629f2..6a50621 --- 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 -- cgit v1.2.1