Version 3.5, from a fresh installationv3.5

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)