Version 3.5, from a fresh installationv3.5

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