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 --- CakeCost.R | 31 +++- CakeErrMin.R | 3 +- CakeFit.R | 205 ++++++++++++++++++++++++ CakeHelpers.R | 188 ++++++++++++++-------- CakeInit.R | 5 +- CakeIrlsFit.R | 394 ++++++++++++++++----------------------------- CakeMcmcFit.R | 482 ++++++++++++++++++++++++++------------------------------ CakeModel.R | 3 +- CakeNullPlot.R | 8 +- CakeOdeSolve.R | 3 +- CakeOlsFit.R | 185 ++++++++-------------- CakePenalties.R | 4 +- CakePlot.R | 3 +- CakeSolutions.R | 12 +- CakeSummary.R | 3 +- 15 files changed, 823 insertions(+), 706 deletions(-) mode change 100644 => 100755 CakeCost.R mode change 100644 => 100755 CakeErrMin.R create mode 100755 CakeFit.R mode change 100644 => 100755 CakeHelpers.R mode change 100644 => 100755 CakeInit.R mode change 100644 => 100755 CakeIrlsFit.R mode change 100644 => 100755 CakeMcmcFit.R mode change 100644 => 100755 CakeModel.R mode change 100644 => 100755 CakeNullPlot.R mode change 100644 => 100755 CakeOdeSolve.R mode change 100644 => 100755 CakeOlsFit.R mode change 100644 => 100755 CakePenalties.R mode change 100644 => 100755 CakePlot.R mode change 100644 => 100755 CakeSolutions.R mode change 100644 => 100755 CakeSummary.R diff --git a/CakeCost.R b/CakeCost.R old mode 100644 new mode 100755 index 6b9fcb3..110275e --- 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 old mode 100644 new mode 100755 index ff49ad5..029db00 --- 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 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 diff --git a/CakeInit.R b/CakeInit.R old mode 100644 new mode 100755 index 57dbc9e..3214e02 --- 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 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 diff --git a/CakeMcmcFit.R b/CakeMcmcFit.R old mode 100644 new mode 100755 index be57cef..0a5bf4a --- 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 old mode 100644 new mode 100755 index 116ef6b..440945c --- 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 old mode 100644 new mode 100755 index 52245ce..f6a6955 --- 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 old mode 100644 new mode 100755 index 263e57c..f1c83f7 --- 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 old mode 100644 new mode 100755 index e9b8f3c..66bf9b4 --- 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 old mode 100644 new mode 100755 index b4f24d4..33d2e0b --- 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 . # 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 old mode 100644 new mode 100755 index da69ea8..63dc76a --- 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 old mode 100644 new mode 100755 index c9e5b88..ecbfa59 --- 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 old mode 100644 new mode 100755 index 5849a4a..21b4ae7 --- 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 -- cgit v1.2.1