diff options
Diffstat (limited to 'CakeCost.R')
| -rwxr-xr-x[-rw-r--r--] | CakeCost.R | 31 |
1 files changed, 23 insertions, 8 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) |