summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rwxr-xr-x[-rw-r--r--]CakeCost.R31
-rwxr-xr-x[-rw-r--r--]CakeErrMin.R3
-rwxr-xr-xCakeFit.R205
-rwxr-xr-x[-rw-r--r--]CakeHelpers.R188
-rwxr-xr-x[-rw-r--r--]CakeInit.R5
-rwxr-xr-x[-rw-r--r--]CakeIrlsFit.R394
-rwxr-xr-x[-rw-r--r--]CakeMcmcFit.R482
-rwxr-xr-x[-rw-r--r--]CakeModel.R3
-rwxr-xr-x[-rw-r--r--]CakeNullPlot.R8
-rwxr-xr-x[-rw-r--r--]CakeOdeSolve.R3
-rwxr-xr-x[-rw-r--r--]CakeOlsFit.R185
-rwxr-xr-x[-rw-r--r--]CakePenalties.R4
-rwxr-xr-x[-rw-r--r--]CakePlot.R3
-rwxr-xr-x[-rw-r--r--]CakeSolutions.R12
-rwxr-xr-x[-rw-r--r--]CakeSummary.R3
15 files changed, 823 insertions, 706 deletions
diff --git a/CakeCost.R b/CakeCost.R
index 6b9fcb3..110275e 100644..100755
--- a/CakeCost.R
+++ b/CakeCost.R
@@ -5,7 +5,8 @@
# Call to approx is only performed if there are multiple non NA values
# which should prevent most of the crashes - Rob Nelson (Tessella)
#
-# 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
@@ -24,6 +25,10 @@
CakeCost <- function (model, obs, x = "time", y = NULL, err = NULL,
weight = "none", scaleVar = FALSE, cost = NULL, ...) {
+ ## Sometimes a fit is encountered for which the model is unable to calculate
+ ## values on the full range of observed values. In this case, we will return
+ ## an infinite cost to ensure this value is not selected.
+ modelCalculatedFully <- all(unlist(obs[x]) %in% unlist(model[x]))
## convert vector to matrix
if (is.vector(obs)) {
@@ -125,7 +130,7 @@ CakeCost <- function (model, obs, x = "time", y = NULL, err = NULL,
for (i in ilist) { # for each observed variable ...
ii <- which(ModNames == Names[i])
if (length(ii) == 0) stop(paste("observed variable not found in model output", Names[i]))
- yMod <- model[,ii]
+ yMod <- model[, ii]
if (type == 2) { # table format
iDat <- which (obs[,1] == Names[i])
if (length(ix) > 0) xDat <- obs[iDat, ix]
@@ -174,17 +179,28 @@ CakeCost <- function (model, obs, x = "time", y = NULL, err = NULL,
if (scaleVar)
Scale <- 1/length(obsdat)
else Scale <- 1
- Res <- (ModVar- obsdat)
- res <- Res/Err
- resScaled <- res*Scale
+ if(!modelCalculatedFully){ # In this case, the model is unable to predict on the full range, set cost to Inf
+ xDat <- 0
+ obsdat <- 0
+ ModVar <- Inf
+ Res <- Inf
+ res <- Inf
+ weight_for_residual <- Inf
+ } else{
+ Res <- (ModVar- obsdat)
+ res <- Res / Err
+ weight_for_residual <- 1 / Err
+ }
+
+ resScaled <- res * Scale
Residual <- rbind(Residual,
data.frame(
name = Names[i],
x = xDat,
obs = obsdat,
mod = ModVar,
- weight = 1/Err,
+ weight = weight_for_residual,
res.unweighted = Res,
res = res))
@@ -201,7 +217,6 @@ CakeCost <- function (model, obs, x = "time", y = NULL, err = NULL,
## SSR
Cost <- sum(CostVar$SSR * CostVar$scale)
-
Lprob <- -sum(log(pmax(0, dnorm(Residual$mod, Residual$obs, Err)))) # avoid log of negative values
#Lprob <- -sum(log(pmax(.Machine$double.xmin, dnorm(Residual$mod, Residual$obs, Err)))) #avoid log(0)
@@ -336,7 +351,7 @@ CakeInternalCostFunctions <- function(mkinmod, state.ini.optim, state.ini.optim.
# HACK to make nls.lm respect the penalty, as it just uses residuals and ignores the cost
if(mC$penalties > 0){
- mC$residuals$res <- mC$residuals$res + mC$penalties / length(mC$residuals$res)
+ mC$residuals$res <- mC$residuals$res + (sign(mC$residuals$res) * mC$penalties / length(mC$residuals$res))
}
return(mC)
diff --git a/CakeErrMin.R b/CakeErrMin.R
index ff49ad5..029db00 100644..100755
--- a/CakeErrMin.R
+++ b/CakeErrMin.R
@@ -3,7 +3,8 @@
## -----------------------------------------------------------------------------
# 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
diff --git a/CakeFit.R b/CakeFit.R
new file mode 100755
index 0000000..7e6c15a
--- /dev/null
+++ b/CakeFit.R
@@ -0,0 +1,205 @@
+# Performs the optimisation routine corresponding to the given optimiser selected.
+# Remark: this function was originally based on the "mkinfit" function, version 0.1.
+#
+# optimiser: A character type indicating which optimiser should be used, one of 'OLS', 'IRLS', 'MCMC'
+# 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).
+# niter: The number of MCMC iterations to apply.
+# atol: The tolerance to apply to the ODE solver.
+# dfopDtMaxIter: The maximum number of iterations to apply to DFOP DT calculation.
+# control: ...
+# useExtraSolver: Whether to use the extra solver for this fit.
+CakeFit <- function(optimiser,
+ cake.model,
+ observed,
+ parms.ini,
+ state.ini,
+ lower = 0,
+ upper = Inf,
+ fixed_parms = NULL,
+ fixed_initials = names(cake.model$diffs)[-1],
+ quiet = FALSE,
+ niter = 1000,
+ verbose = TRUE,
+ seed = NULL,
+ atol = 1e-6,
+ dfopDtMaxIter = 10000,
+ control = list(),
+ useExtraSolver = FALSE) {
+ env = environment()
+
+ # Get model variable names
+ mod_vars <- names(cake.model$diffs)
+
+ # Subset dataframe with mapped (modelled) variables
+ observed <- subset(observed, name %in% names(cake.model$map))
+
+ ERR <- rep(1, nrow(observed))
+ observed <- cbind(observed, err = ERR)
+
+ # Get names of observed variables
+ obs_vars <- unique(as.character(observed$name))
+
+ # 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
+ # 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) {
+ # numerical method, replace with analytical for parent-only fits
+ 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)))
+ }
+ }
+
+ fitted_with_extra_solver <- FALSE
+
+ optimiserSpecificSetup <- GetOptimiserSpecificSetup(optimiser)
+ optimiserSpecificSetupOutput <- optimiserSpecificSetup(cake.model,
+ state.ini.optim,
+ state.ini.optim.boxnames,
+ state.ini.fixed,
+ parms.fixed,
+ observed,
+ mkindiff = mkindiff,
+ quiet,
+ atol = atol,
+ solution = solution,
+ control = control,
+ seed = seed)
+ list2env(optimiserSpecificSetupOutput, env)
+
+ pnames <- names(c(state.ini.optim, parms.optim))
+ costForExtraSolver <- function(P) {
+ names(P) <- pnames
+ FF <<- costFunctions$cost(P)
+ return(FF$model)
+ }
+
+ optimiserFitRoutine <- GetOptimisationRoutine(optimiser)
+ optimiserFitResult <- optimiserFitRoutine(costFunctions,
+ costForExtraSolver,
+ useExtraSolver,
+ c(state.ini.optim, parms.optim),
+ lower,
+ upper,
+ control,
+ cake.model = cake.model,
+ tol = tol,
+ fitted_with_extra_solver = fitted_with_extra_solver,
+ observed = observed,
+ maxIter = maxIter,
+ costWithStatus = costWithStatus,
+ niter = niter,
+ verbose = verbose)
+ list2env(optimiserFitResult, env)
+
+ optimiserSpecificWrapUp <- GetOptimiserSpecificWrapUp(optimiser)
+ optimiserSpecificWrapUpOutput <- optimiserSpecificWrapUp(fit,
+ parms.fixed = parms.fixed,
+ state.ini = state.ini,
+ observed = observed,
+ res = res,
+ seed = seed,
+ costWithStatus = costWithStatus)
+ list2env(optimiserSpecificWrapUpOutput, env)
+
+ # We need to return some more data for summary and plotting
+ fit$solution <- solution
+
+ if (solution == "deSolve") {
+ 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 = solution, atol)
+ out_transformed <- PostProcessOdeOutput(out_predicted, cake.model, atol)
+
+ predicted_long <- wide_to_long(out_transformed, time = "time")
+ fit$predicted <- out_transformed
+
+ # 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)
+
+ # Calculate dissipation times DT50 and DT90 and formation fractions
+ 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)
+ }
+
+ fit$ioreRepDT <- CakeIORERepresentativeDT("Parent", parms.all)
+ fit$fomcRepDT <- CakeFOMCBackCalculatedDT(parms.all)
+ fit$penalties <- CakePenaltiesLong(parms.all, out_transformed, observed)
+
+ # Collect observed, predicted and residuals
+ 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$costFunctions <- costFunctions
+ fit$upper <- upper
+ fit$lower <- lower
+
+ fit$atol <- atol
+
+ if (optimiser == "MCMC") {
+ sq <- data$residual * data$residual
+ fit$ssr <- sum(sq)
+ }
+
+ return(fit)
+} \ No newline at end of file
diff --git a/CakeHelpers.R b/CakeHelpers.R
index cc69509..0e9c518 100644..100755
--- 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
diff --git a/CakeInit.R b/CakeInit.R
index 57dbc9e..3214e02 100644..100755
--- a/CakeInit.R
+++ b/CakeInit.R
@@ -4,7 +4,7 @@
# Call with chdir=TRUE so it can load our source from the current directory
# Tessella Project Reference: 6245, 7247, 8361, 7414, 10091
-# Copyright (C) 2011-2020 Syngenta
+# Copyright (C) 2011-2022 Syngenta
# 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
@@ -34,13 +34,14 @@ source("CakeErrMin.R")
source("CakeHelpers.R")
source("CakeSolutions.R")
source("CakeOdeSolve.R")
+source("CakeFit.R")
options(width=1000)
options(error=recover)
options(warn=-1)
# When running from C#, this must match the C# CAKE version
-CAKE.version<-"3.4"
+CAKE.version<-"3.5"
# The number of data points to use to draw lines on CAKE plots.
CAKE.plots.resolution<-401
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
diff --git a/CakeMcmcFit.R b/CakeMcmcFit.R
index be57cef..0a5bf4a 100644..100755
--- a/CakeMcmcFit.R
+++ b/CakeMcmcFit.R
@@ -1,6 +1,7 @@
# Some of the CAKE R modules are based on mkin,
# Based on mcmckinfit as modified by Bayer
-# 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
@@ -32,270 +33,241 @@
# dfopDtMaxIter: The maximum number of iterations to apply to DFOP DT calculation.
# control: ...
# useExtraSolver: Whether to use the extra solver for this fit (only used for the initial first fit).
-CakeMcmcFit <- function (cake.model,
- observed,
- parms.ini,
- state.ini,
- lower,
- upper,
- fixed_parms = NULL,
- fixed_initials,
- quiet = FALSE,
- niter = 1000,
- verbose = TRUE,
- seed=NULL,
- 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)
- 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]])))
+CakeMcmcFit <- function(cake.model,
+ observed,
+ parms.ini,
+ state.ini,
+ lower,
+ upper,
+ fixed_parms = NULL,
+ fixed_initials,
+ quiet = FALSE,
+ niter = 1000,
+ verbose = TRUE,
+ seed = NULL,
+ atol = 1e-6,
+ dfopDtMaxIter = 10000,
+ control = list(),
+ useExtraSolver = FALSE) {
+
+ fit <- CakeFit("MCMC",
+ cake.model,
+ observed,
+ parms.ini,
+ state.ini,
+ lower,
+ upper,
+ fixed_parms,
+ fixed_initials,
+ quiet,
+ niter = niter,
+ verbose = verbose,
+ seed = seed,
+ atol = atol,
+ dfopDtMaxIter = dfopDtMaxIter,
+ control = control,
+ useExtraSolver = useExtraSolver)
+
+ return(fit)
+}
+
+GetMcmcSpecificSetup <- function() {
+ return(function(
+ cake.model,
+ state.ini.optim,
+ state.ini.optim.boxnames,
+ state.ini.fixed,
+ parms.fixed,
+ observed,
+ mkindiff,
+ quiet,
+ atol,
+ solution,
+ ...) {
+ seed <- list(...)$seed
+
+ costFunctions <- CakeInternalCostFunctions(cake.model, state.ini.optim, state.ini.optim.boxnames, state.ini.fixed,
+ parms.fixed, observed, mkindiff = mkindiff, quiet, atol = atol, solution = solution)
+
+ bestIteration <<- -1;
+ costWithStatus <- function(P, ...) {
+ r <- costFunctions$cost(P)
+ if (r$cost == costFunctions$get.best.cost()) {
+ bestIteration <<- costFunctions$get.calls();
+ cat(' MCMC best so far: c', r$cost, 'it:', bestIteration, '\n')
+ }
+
+ arguments <- list(...)
+ if (costFunctions$get.calls() <= arguments$maxCallNo) {
+ cat("MCMC Call no.", costFunctions$get.calls(), '\n')
+ }
+
+ return(r)
}
-
- 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)
- bestIteration <- -1;
- costWithStatus <- function(P, ...){
- r <- costFunctions$cost(P)
- if(r$cost == costFunctions$get.best.cost()) {
- bestIteration<<-costFunctions$get.calls();
- cat(' MCMC best so far: c', r$cost, 'it:', bestIteration, '\n')
+
+ # Set the seed
+ if (is.null(seed)) {
+ # No seed so create a random one so there is something to report
+ seed <- runif(1, 0, 2 ^ 31 - 1)
}
-
- arguments <- list(...)
- if (costFunctions$get.calls() <= arguments$maxCallNo)
- {
- cat("MCMC Call no.", costFunctions$get.calls(), '\n')
+
+ seed <- as.integer(seed)
+ set.seed(seed)
+
+ return(list(costFunctions = costFunctions, costWithStatus = costWithStatus, maxIter = NULL, tol = NULL, seed = seed))
+ })
+}
+
+GetMcmcOptimisationRoutine <- function() {
+ return(function(costFunctions, costForExtraSolver, useExtraSolver, parms, lower, upper, control, ...) {
+ mcmcArgs <- list(...)
+ cake.model <- mcmcArgs$cake.model
+ costWithStatus <- mcmcArgs$costWithStatus
+ observed <- mcmcArgs$observed
+ niter <- mcmcArgs$niter
+ verbose <- mcmcArgs$verbose
+
+ # Runs a pre-fit with no weights first, followed by a weighted step (extra solver with first, not with second)
+
+ # Run optimiser, no weighting
+ fitStepResult <- RunFitStep(costFunctions$cost, costForExtraSolver, useExtraSolver, parms, lower, upper, control)
+ fit <- fitStepResult$fit
+ fitted_with_extra_solver <- fitStepResult$fitted_with_extra_solver
+
+
+ # Process extra solver output if it was used
+ if (fitted_with_extra_solver) {
+ fit <- GetFitValuesAfterExtraSolver(fit, FF)
}
- return(r)
- }
+ # One reweighted estimation
+ # Estimate the error variance(sd)
+ tmpres <- fit$residuals
+ oldERR <- observed$err
+ err <- rep(NA, length(cake.model$map))
- # Set the seed
- if ( is.null(seed) ) {
- # No seed so create a random one so there is something to report
- seed <- runif(1,0,2^31-1)
- }
-
- seed <- as.integer(seed)
- set.seed(seed)
-
- ## ############# Get Initial Paramtervalues #############
- ## Start with no weighting
- if(useExtraSolver)
- {
- a <- try(fit <- modFit(costFunctions$cost, c(state.ini.optim, parms.optim), lower = lower, upper = upper, method='Port',control=control))
-
- 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)
+ for (i in 1:length(cake.model$map)) {
+ box <- names(cake.model$map)[i]
+ ind <- which(names(tmpres) == box)
+ tmp <- tmpres[ind]
+ err[i] <- sd(tmp)
}
- }
- 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,...)
- }
-
- ## ############## One reweighted estimation ############################
- ## Estimate the error variance(sd)
- tmpres <- fit$residuals
- oldERR <- observed$err
- err <- rep(NA,length(cake.model$map))
-
- for(i in 1:length(cake.model$map)) {
- box <- names(cake.model$map)[i]
- ind <- which(names(tmpres)==box)
- tmp <- tmpres[ind]
- err[i] <- sd(tmp)
- }
-
- names(err) <- names(cake.model$map)
- ERR <- err[as.character(observed$name)]
- observed$err <-ERR
- costFunctions$set.error(ERR)
-
- olderr <- rep(1,length(cake.model$map))
- diffsigma <- sum((err-olderr)^2)
- ## At least do one iteration step to get a weighted LS
- fit <- modFit(costFunctions$cost, fit$par, lower = lower, upper = upper, ...)
- ## Use this as the Input for MCMC algorithm
- ## ##########################
- fs <- summary(fit)
- cov0 <- if(all(is.na(fs$cov.scaled))) NULL else fs$cov.scaled*2.4^2/length(fit$par)
- var0 <- fit$var_ms_unweighted
- costFunctions$set.calls(0); costFunctions$reset.best.cost()
- res <- modMCMC(costWithStatus, maxCallNo=niter, fit$par,...,jump=cov0,lower=lower,upper=upper,prior=NULL,var0=var0,wvar0=0.1,niter=niter,outputlength=niter,burninlength=0,updatecov=niter,ntrydr=1,drscale=NULL,verbose=verbose)
-
- # Construct the fit object to return
- 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)))
-
- fit$mkindiff <- mkindiff
- fit$map <- cake.model$map
- fit$diffs <- cake.model$diffs
-
- # Replace mean from modFit with mean from modMCMC
- fnm <- function(x) mean(res$pars[,x])
- fit$par <- sapply(dimnames(res$pars)[[2]],fnm)
- fit$bestpar <- res$bestpar
- fit$costfn <- costWithStatus
-
- # Disappearence times
- parms.all <- c(fit$par, parms.fixed)
- obs_vars <- unique(as.character(observed$name))
- 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)
- }
-
- fit$ioreRepDT <- CakeIORERepresentativeDT("Parent", parms.all)
- fit$fomcRepDT <- CakeFOMCBackCalculatedDT(parms.all)
-
- # Ensure initial state is at time 0
- obstimes <- unique(c(0,observed$time))
-
- # Solve the system
- out_predicted <- CakeOdeSolve(fit, obstimes, solution = "deSolve", atol)
- out_transformed <- PostProcessOdeOutput(out_predicted, cake.model, atol)
-
- fit$predicted <- out_transformed
- fit$penalties <- CakePenaltiesLong(parms.all, out_transformed, observed)
-
- predicted_long <- wide_to_long(out_transformed,time="time")
- obs_vars <- unique(as.character(observed$name))
- fit$errmin <- CakeChi2(cake.model, observed, predicted_long, obs_vars, parms.optim, state.ini.optim, state.ini, parms.ini, fit$fixed)
-
- 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
-
- sq <- data$residual * data$residual
- fit$ssr <- sum(sq)
-
- fit$seed <- seed
-
- fit$res <- res
- class(fit) <- c("CakeMcmcFit", "mkinfit", "modFit")
- return(fit)
+
+ names(err) <- names(cake.model$map)
+ ERR <- err[as.character(observed$name)]
+ observed$err <- ERR
+ costFunctions$set.error(ERR)
+
+ olderr <- rep(1, length(cake.model$map))
+ diffsigma <- sum((err - olderr) ^ 2)
+
+ ## At least do one iteration step to get a weighted LS
+ fit <- modFit(f = costFunctions$cost, p = fit$par, lower = lower, upper = upper)
+
+ # Run MCMC optimiser with output from weighted fit
+ # Apply iterative re-weighting here (iterations fixed to 1 for now):
+ # Do modMCMC as below, we should also pass in the final priors for subsequent iterations
+ # Use modMCMC average as input to next modMCMC run (as done in block 3 to get final parameters)
+ # use errors from previous step as inputs to modMCMC cov0 and var0 at each iteration
+
+ fs <- summary(fit)
+ cov0 <- if (all(is.na(fs$cov.scaled))) NULL else fs$cov.scaled * 2.4 ^ 2 / length(fit$par)
+ var0 <- fit$var_ms_unweighted
+ costFunctions$set.calls(0);
+ costFunctions$reset.best.cost()
+ res <- modMCMC(f = costWithStatus, p = fit$par, maxCallNo = niter, jump = cov0, lower = lower, upper = upper, prior = NULL, var0 = var0, wvar0 = 0.1, niter = niter, outputlength = niter, burninlength = 0, updatecov = niter, ntrydr = 1, drscale = NULL, verbose = verbose)
+ return(list(fit = fitStepResult$fit, fitted_with_extra_solver = fitStepResult$fitted_with_extra_solver, res = res))
+ })
}
+GetMcmcSpecificWrapUp <- function() {
+ return(function(fit, ...) {
+ args <- list(...)
+ res <- args$res
+ seed <- args$seed
+ costWithStatus <- args$costWithStatus
+ observed <- args$observed
+ parms.fixed <- args$parms.fixed
+
+ # Replace mean from modFit with mean from modMCMC
+ fnm <- function(x) mean(res$pars[, x])
+ fit$par <- sapply(dimnames(res$pars)[[2]], fnm)
+ fit$bestpar <- res$bestpar
+ fit$costfn <- costWithStatus
+ parms.all <- c(fit$par, parms.fixed)
+
+ data <- observed
+ data$err <- rep(NA, length(data$time))
+
+ fit$seed <- seed
+ fit$res <- res
+
+ np <- length(parms.all)
+ fit$rank <- np
+ fit$df.residual <- length(fit$residuals) - fit$rank
+
+ class(fit) <- c("CakeMcmcFit", "mkinfit", "modFit") # Note different class to other optimisers
+ return(list(fit = fit, parms.all = parms.all, data = data))
+ })
+}
# Summarise a fit
-summary.CakeMcmcFit <- function(object, data = TRUE, distimes = TRUE, halflives = TRUE, ff = TRUE, cov = FALSE,...) {
- param <- object$par
- pnames <- names(param)
- p <- length(param)
- #covar <- try(solve(0.5*object$hessian), silent = TRUE) # unscaled covariance
- mcmc <- object$res
- covar <- cov(mcmc$pars)
-
- rdf <- object$df.residual
-
+# The MCMC summary is separate from the others due to the difference in the outputs of the modMCMC and modFit.
+summary.CakeMcmcFit <- function(object, data = TRUE, distimes = TRUE, halflives = TRUE, ff = TRUE, cov = FALSE, ...) {
+ param <- object$par
+ pnames <- names(param)
+ p <- length(param)
+ #covar <- try(solve(0.5*object$hessian), silent = TRUE) # unscaled covariance
+ mcmc <- object$res
+ covar <- cov(mcmc$pars)
+
+ rdf <- object$df.residual
+
message <- "ok"
- rownames(covar) <- colnames(covar) <-pnames
-
+ rownames(covar) <- colnames(covar) <- pnames
+
#se <- sqrt(diag(covar) * resvar)
- fnse <- function(x) sd(mcmc$pars[,x]) #/sqrt(length(mcmc$pars[,x]))
- se <- sapply(dimnames(mcmc$pars)[[2]],fnse)
-
- tval <- param / se
+ fnse <- function(x) sd(mcmc$pars[, x]) #/sqrt(length(mcmc$pars[,x]))
+ se <- sapply(dimnames(mcmc$pars)[[2]], fnse)
+
+ tval <- param / se
- if (!all(object$start$lower >=0)) {
- message <- "Note that the one-sided t-test may not be appropriate if
- parameter values below zero are possible."
- warning(message)
- } else message <- "ok"
+ if (!all(object$start$lower >= 0)) {
+ message <- "Note that the one-sided t-test may not be appropriate if
+ parameter values below zero are possible."
+ warning(message)
+ } else message <- "ok"
# Filter the values for t-test, only apply t-test to k-values
- t.names <- grep("k(\\d+)|k_(.*)", names(tval), value = TRUE)
- t.rest <- setdiff(names(tval), t.names)
+ t.names <- grep("k(\\d+)|k_(.*)", names(tval), value = TRUE)
+ t.rest <- setdiff(names(tval), t.names)
t.values <- c(tval)
t.values[t.rest] <- NA
t.result <- pt(t.values, rdf, lower.tail = FALSE)
-
+
# Now set the values we're not interested in for the lower
# and upper bound we're not interested in to NA
t.param <- c(param)
t.param[t.names] <- NA
# calculate the 90% confidence interval
alpha <- 0.10
- lci90 <- t.param + qt(alpha/2,rdf)*se
- uci90 <- t.param + qt(1-alpha/2,rdf)*se
-
+ lci90 <- t.param + qt(alpha / 2, rdf) * se
+ uci90 <- t.param + qt(1 - alpha / 2, rdf) * se
+
# calculate the 95% confidence interval
alpha <- 0.05
- lci95 <- t.param + qt(alpha/2,rdf)*se
- uci95 <- t.param + qt(1-alpha/2,rdf)*se
+ lci95 <- t.param + qt(alpha / 2, rdf) * se
+ uci95 <- t.param + qt(1 - alpha / 2, rdf) * se
param <- cbind(param, se, tval, t.result, lci90, uci90, lci95, uci95)
dimnames(param) <- list(pnames, c("Estimate", "Std. Error",
"t value", "Pr(>t)", "Lower CI (90%)", "Upper CI (90%)", "Lower CI (95%)", "Upper CI (95%)"))
-
- # Residuals from mean of MCMC fit
- resvar <- object$ssr/ rdf
- modVariance <- object$ssr / length(object$data$residual)
-
- ans <- list(ssr = object$ssr,
+
+ # Residuals from mean of MCMC fit
+ resvar <- object$ssr / rdf
+ modVariance <- object$ssr / length(object$data$residual)
+
+ ans <- list(ssr = object$ssr,
residuals = object$data$residuals,
residualVariance = resvar,
sigma = sqrt(resvar),
@@ -305,24 +277,24 @@ summary.CakeMcmcFit <- function(object, data = TRUE, distimes = TRUE, halflives
info = object$info, niter = object$iterations,
stopmess = message,
par = param)
-
- ans$diffs <- object$diffs
- ans$data <- object$data
- ans$additionalstats <- CakeAdditionalStats(object$data)
- ans$seed <- object$seed
-
- ans$start <- object$start
- ans$fixed <- object$fixed
- ans$errmin <- object$errmin
- ans$penalties <- object$penalties
- if(distimes){
- ans$distimes <- object$distimes
- ans$extraDT50 <- object$extraDT50
- ans$ioreRepDT <- object$ioreRepDT
- ans$fomcRepDT <- object$fomcRepDT
- }
- if(halflives) ans$halflives <- object$halflives
- if(ff) ans$ff <- object$ff
- class(ans) <- c("summary.CakeFit","summary.mkinfit", "summary.modFit")
- return(ans)
+
+ ans$diffs <- object$diffs
+ ans$data <- object$data
+ ans$additionalstats <- CakeAdditionalStats(object$data)
+ ans$seed <- object$seed
+
+ ans$start <- object$start
+ ans$fixed <- object$fixed
+ ans$errmin <- object$errmin
+ ans$penalties <- object$penalties
+ if (distimes) {
+ ans$distimes <- object$distimes
+ ans$extraDT50 <- object$extraDT50
+ ans$ioreRepDT <- object$ioreRepDT
+ ans$fomcRepDT <- object$fomcRepDT
+ }
+ if (halflives) ans$halflives <- object$halflives
+ if (ff) ans$ff <- object$ff
+ class(ans) <- c("summary.CakeFit", "summary.mkinfit", "summary.modFit")
+ return(ans)
}
diff --git a/CakeModel.R b/CakeModel.R
index 116ef6b..440945c 100644..100755
--- a/CakeModel.R
+++ b/CakeModel.R
@@ -4,7 +4,8 @@
# Portions Johannes Ranke, 2010
# Contact: mkin-devel@lists.berlios.de
-# 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
diff --git a/CakeNullPlot.R b/CakeNullPlot.R
index 52245ce..f6a6955 100644..100755
--- a/CakeNullPlot.R
+++ b/CakeNullPlot.R
@@ -2,7 +2,8 @@
# Generates model curves so the GUI can plot them. Used when all params are fixed so there was no fit
# Some of the CAKE R modules are based on mkin
# Based on code in IRLSkinfit
-# 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
@@ -118,6 +119,11 @@ CakeNullPlot <- function(model, parms.ini, state.ini, observed, xlim = range(obs
data<-data[order(data$variable,data$time),]
cat("\nData:\n")
print(format(data, digits = digits, scientific = FALSE,...), row.names = FALSE)
+
+ additionalStats <- CakeAdditionalStats(data)
+
+ cat("\nAdditional Stats:\n");
+ print(additionalStats, digits = digits, ...)
# Solve at output times
out <- ode(
diff --git a/CakeOdeSolve.R b/CakeOdeSolve.R
index 263e57c..f1c83f7 100644..100755
--- a/CakeOdeSolve.R
+++ b/CakeOdeSolve.R
@@ -3,7 +3,8 @@
## -----------------------------------------------------------------------------
# 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
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
diff --git a/CakePenalties.R b/CakePenalties.R
index b4f24d4..33d2e0b 100644..100755
--- a/CakePenalties.R
+++ b/CakePenalties.R
@@ -1,6 +1,7 @@
# $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
@@ -17,6 +18,7 @@
# along with this program. If not, see <http://www.gnu.org/licenses/>.
# Penalty function for fits that should not be accepted
+# Note: If the penalty scale is updated from 10 then the reported % penalty in ReportGenerator.cs will need to be updated
CakePenalties<-function(params, modelled, obs, penalty.scale = 10, ...){
# Because the model cost is roughly proportional to points,
# make the penalties scale that way too for consistent behaviour
diff --git a/CakePlot.R b/CakePlot.R
index da69ea8..63dc76a 100644..100755
--- a/CakePlot.R
+++ b/CakePlot.R
@@ -1,7 +1,8 @@
#$Id$
# Generates fitted curves so the GUI can plot them
# Based on code in IRLSkinfit
-# 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
diff --git a/CakeSolutions.R b/CakeSolutions.R
index c9e5b88..ecbfa59 100644..100755
--- a/CakeSolutions.R
+++ b/CakeSolutions.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
+# 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
@@ -38,5 +39,12 @@ HS.solution <- function (t, parent.0, k1, k2, tb) {
# Produces solutions to the IORE equation given times t and parameters M_0, k and N.
IORE.solution <- function(t, parent.0, k, N) {
- parent = (parent.0^(1 - N) - (1 - N) * k * t)^(1/(1-N))
+ if (N == 1) {
+ parent = parent.0 * exp(-k * t)
+ }
+ else {
+ parent = (parent.0 ^ (1 - N) - (1 - N) * k * t) ^ (1 / (1 - N))
+ }
+
+ return(parent)
} \ No newline at end of file
diff --git a/CakeSummary.R b/CakeSummary.R
index 5849a4a..21b4ae7 100644..100755
--- a/CakeSummary.R
+++ b/CakeSummary.R
@@ -3,7 +3,8 @@
# and display the summary itself.
# Parts modified from package mkin
-# Parts Tessella for Syngenta: Copyright (C) 2011-2020 Syngenta
+# Parts 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

Contact - Imprint