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 --- CakeHelpers.R | 188 +++++++++++++++++++++++++++++++++++++++------------------- 1 file changed, 126 insertions(+), 62 deletions(-) mode change 100644 => 100755 CakeHelpers.R (limited to 'CakeHelpers.R') diff --git a/CakeHelpers.R b/CakeHelpers.R old mode 100644 new mode 100755 index cc69509..0e9c518 --- a/CakeHelpers.R +++ b/CakeHelpers.R @@ -1,6 +1,7 @@ # $Id$ # Some of the CAKE R modules are based on mkin, -# Developed by Tessella Ltd for Syngenta: Copyright (C) 2011-2020 Syngenta +# 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 @@ -19,13 +20,13 @@ # 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 + 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 + parametersOnUpperBound = which(parameters == upper) + parameters[parametersOnUpperBound] <- parameters[parametersOnUpperBound] * (1 - .Machine$double.neg.eps) - .Machine$double.xmin - return(parameters) + return(parameters) } # Adjusts stated initial values to put into the ODE solver. @@ -36,25 +37,25 @@ ShiftAwayFromBoundaries <- function(parameters, lower, upper) { # # 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]]) - } + 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)]) + + # 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. @@ -65,34 +66,34 @@ AdjustOdeInitialValues <- function(odeini, cake.model, odeparms) { # # 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 + 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 } - - # 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]]]) + + # 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) + + return(out_transformed) } # Reorganises data in a wide format to a long format. @@ -102,17 +103,80 @@ PostProcessOdeOutput <- function(odeoutput, cake.model, atol) { # # 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])), + 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) + + 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()) } \ No newline at end of file -- cgit v1.2.1