Version 3.5, from a fresh installationv3.5

1 files changed, 126 insertions, 62 deletions

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