Version 3.4 as just publicly announcedv3.4

Peter Rainbird just announced the release on the PFMODELS email list.

14 files changed, 442 insertions, 321 deletions

diff --git a/CakeCost.R b/CakeCost.R
index 8fb94ef..6b9fcb3 100644
--- a/CakeCost.R
+++ b/CakeCost.R
@@ -5,8 +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-2016 Syngenta
-# Tessella Project Reference: 6245, 7247, 8361, 7414
+# Modifications developed by Tessella for Syngenta, Copyright (C) 2011-2020 Syngenta
+# Tessella Project Reference: 6245, 7247, 8361, 7414, 10091
 #
 # The CAKE R modules are free software: you can
 # redistribute them and/or modify them under the
diff --git a/CakeErrMin.R b/CakeErrMin.R
index 9a074a4..ff49ad5 100755..100644
--- a/CakeErrMin.R
+++ b/CakeErrMin.R
@@ -3,8 +3,8 @@
 ## -----------------------------------------------------------------------------
 # Some of the CAKE R modules are based on mkin.
 #
-# Modifications developed by Tessella for Syngenta, Copyright (C) 2011-2016 Syngenta
-# Tessella Project Reference: 6245, 7247, 8361, 7414
+# Modifications developed by Tessella for Syngenta, Copyright (C) 2011-2020 Syngenta
+# Tessella Project Reference: 6245, 7247, 8361, 7414, 10091
 #
 # The CAKE R modules are free software: you can
 # redistribute them and/or modify them under the
diff --git a/CakeHelpers.R b/CakeHelpers.R
index 5af51cc..cc69509 100755..100644
--- a/CakeHelpers.R
+++ b/CakeHelpers.R
@@ -1,7 +1,7 @@
 # $Id$
 # Some of the CAKE R modules are based on mkin, 
-# Developed by Tessella Ltd for Syngenta: Copyright (C) 2011-2016 Syngenta
-# Tessella Project Reference: 6245, 7247, 8361, 7414
+# Developed by Tessella Ltd for Syngenta: Copyright (C) 2011-2020 Syngenta
+# Tessella Project Reference: 6245, 7247, 8361, 7414, 10091
 
 #    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
diff --git a/CakeInit.R b/CakeInit.R
index 72cb433..57dbc9e 100644
--- a/CakeInit.R
+++ b/CakeInit.R
@@ -2,9 +2,9 @@
 
 # Purpose: Load the modules required by CAKE
 # Call with chdir=TRUE so it can load our source from the current directory
-# Tessella Project Reference: 6245, 7247, 8361, 7414
+# Tessella Project Reference: 6245, 7247, 8361, 7414, 10091
 
-# Copyright (C) 2011-2016 Syngenta
+# Copyright (C) 2011-2020 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
@@ -33,13 +33,14 @@ source("CakeSummary.R")
 source("CakeErrMin.R")
 source("CakeHelpers.R")
 source("CakeSolutions.R")
+source("CakeOdeSolve.R")
 
 options(width=1000)
 options(error=recover)
 options(warn=-1)
 
 # When running from C#, this must match the C# CAKE version
-CAKE.version<-"3.3"
+CAKE.version<-"3.4"
 
 # 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 92c45da..f1629f2 100644
--- a/CakeIrlsFit.R
+++ b/CakeIrlsFit.R
@@ -1,8 +1,8 @@
 #$Id$
 #
 # Some of the CAKE R modules are based on mkin
-# Modifications developed by Tessella for Syngenta: Copyright (C) 2011-2016 Syngenta
-# Tessella Project Reference: 6245, 7247, 8361, 7414
+# Modifications developed by Tessella for Syngenta: Copyright (C) 2011-2020 Syngenta
+# Tessella Project Reference: 6245, 7247, 8361, 7414, 10091
  
 #    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
@@ -31,7 +31,7 @@
 # fixed_initials: A vector of compartments with fixed initial concentrations.
 # quiet: Whether the internal cost functions should execute more quietly than normal (less output).
 # atol: The tolerance to apply to the ODE solver.
-# sannMaxIter: The maximum number of iterations to apply to SANN processes.
+# 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, 
@@ -44,7 +44,7 @@ CakeIrlsFit <- function (cake.model,
                          fixed_initials = names(cake.model$diffs)[-1], 
                          quiet = FALSE, 
                          atol=1e-6, 
-                         sannMaxIter = 10000, 
+                         dfopDtMaxIter = 10000, 
                          control=list(),
                          useExtraSolver = FALSE,
                          ...) 
@@ -56,7 +56,7 @@ CakeIrlsFit <- function (cake.model,
     observed <- cbind(observed,err=ERR)
     fitted_with_extra_solver <- 0
     
-    obs_vars = unique(as.character(observed$name))
+    obs_vars <- unique(as.character(observed$name))
     
     if (is.null(names(parms.ini))) {
         names(parms.ini) <- cake.model$parms
@@ -242,45 +242,47 @@ CakeIrlsFit <- function (cake.model,
         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
-        fitForHessian <- modFit(costFunctions$cost, fit$par, lower=lower, upper=upper, method='L-BFGS-B', control=list())
-        
-        fit$solnpHessian <- fit$hessian
-        fit$hessian <- fitForHessian$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
     
-    out_predicted <- costFunctions$get.predicted()
-    
-    predicted_long <- wide_to_long(out_predicted, time = "time")
-    fit$predicted <- out_predicted
     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$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$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_predicted, observed)
+    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,sannMaxIter)
-      fit$extraDT50[compartment.name, ] = CakeExtraDT(type, compartment.name, parms.all)
+      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$ioreRepDT <- CakeIORERepresentativeDT("Parent", parms.all)
+    fit$fomcRepDT <- CakeFOMCBackCalculatedDT(parms.all)
     
     data <- merge(observed, predicted_long, by = c("time", "name"))
     
diff --git a/CakeMcmcFit.R b/CakeMcmcFit.R
index 65c9431..be57cef 100644
--- a/CakeMcmcFit.R
+++ b/CakeMcmcFit.R
@@ -1,7 +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-2016 Syngenta
-# Tessella Project Reference: 6245, 7247, 8361, 7414
+# Modifications developed by Tessella for Syngenta: Copyright (C) 2011-2020 Syngenta
+# Tessella Project Reference: 6245, 7247, 8361, 7414, 10091
 
 #    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
@@ -29,7 +29,7 @@
 # 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.
-# sannMaxIter: The maximum number of iterations to apply to SANN processes.
+# 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, 
@@ -45,7 +45,7 @@ CakeMcmcFit <- function (cake.model,
                          verbose = TRUE, 
                          seed=NULL, 
                          atol=1e-6, 
-                         sannMaxIter = 10000, 
+                         dfopDtMaxIter = 10000, 
                          control=list(),
                          useExtraSolver = FALSE,
                          ...) 
@@ -55,7 +55,7 @@ CakeMcmcFit <- function (cake.model,
     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))
+    obs_vars <- unique(as.character(observed$name))
     
     if (is.null(names(parms.ini)))  {
         names(parms.ini) <- cake.model$parms
@@ -100,17 +100,23 @@ CakeMcmcFit <- function (cake.model,
             bestIteration<<-costFunctions$get.calls();
              cat(' MCMC best so far: c', r$cost, 'it:', bestIteration, '\n')
         }
-        cat("MCMC Call no.", costFunctions$get.calls(), '\n')
+        
+        arguments <- list(...)
+        if (costFunctions$get.calls() <= arguments$maxCallNo)
+        {
+          cat("MCMC Call no.", costFunctions$get.calls(), '\n')
+        }
+
         return(r)
     }
 
     # 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 <- runif(1,0,2^31-1)
     }
     
-    seed = as.integer(seed)
+    seed <- as.integer(seed)
     set.seed(seed)
     
     ## ############# Get Initial Paramtervalues   #############
@@ -161,16 +167,16 @@ CakeMcmcFit <- function (cake.model,
     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,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)
+    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)), 
+    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)), 
+    fit$fixed$type <- c(rep("state", length(state.ini.fixed)), 
         rep("deparm", length(parms.fixed)))
 
     fit$mkindiff <- mkindiff
@@ -185,46 +191,32 @@ CakeMcmcFit <- function (cake.model,
    
     # Disappearence times
     parms.all <- c(fit$par, parms.fixed)
-    obs_vars = unique(as.character(observed$name))
+    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,sannMaxIter)
-        fit$extraDT50[compartment.name, ] = CakeExtraDT(type, compartment.name, parms.all)
+        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)
-    
-    # Fits to data
-    odeini <- state.ini
-    odeparms <- parms.all
-    # If this step modelled the start amount, include it. Lookup in vectors returns NA not null if missing
-    if(!is.na(odeparms["Parent_0"])) { odeini['Parent'] = odeparms["Parent_0"] }
-   
-    # Solve the system
-    evalparse <- function(string)
-    {
-       eval(parse(text=string), as.list(c(odeparms, odeini)))
-    }
+    fit$ioreRepDT <- CakeIORERepresentativeDT("Parent", parms.all)
+    fit$fomcRepDT <- CakeFOMCBackCalculatedDT(parms.all)
   
     # Ensure initial state is at time 0
     obstimes <- unique(c(0,observed$time))
-
-    odeini <- AdjustOdeInitialValues(odeini, cake.model, odeparms)
-   
-    out <- ode(y = odeini, times = obstimes, func = mkindiff, parms = odeparms, atol = atol)
    
-    out_transformed <- PostProcessOdeOutput(out, cake.model, atol)
+    # 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(odeparms, out_transformed, observed)
+    fit$penalties <- CakePenaltiesLong(parms.all, out_transformed, observed)
 
     predicted_long <- wide_to_long(out_transformed,time="time")
-    obs_vars = unique(as.character(observed$name))
+    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
@@ -235,11 +227,12 @@ CakeMcmcFit <- function (cake.model,
     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$seed <- seed
    
     fit$res <- res
     class(fit) <- c("CakeMcmcFit", "mkinfit", "modFit")
@@ -282,7 +275,7 @@ summary.CakeMcmcFit <- function(object, data = TRUE, distimes = TRUE, halflives
     
     # 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 <- c(param)
     t.param[t.names] <- NA
     # calculate the 90% confidence interval
     alpha <- 0.10
@@ -302,24 +295,26 @@ summary.CakeMcmcFit <- function(object, data = TRUE, distimes = TRUE, halflives
    resvar <- object$ssr/ rdf
    modVariance <- object$ssr / length(object$data$residual)
   
-   ans <- list(residuals = object$data$residuals,
-                residualVariance = resvar,
-                sigma = sqrt(resvar),
-                modVariance = modVariance,
-                df = c(p, rdf), cov.unscaled = covar,
-                cov.scaled = covar * resvar,
-                info = object$info, niter = object$iterations,
-                stopmess = message,
-                par = param)
+   ans <- list(ssr = object$ssr,
+               residuals = object$data$residuals,
+               residualVariance = resvar,
+               sigma = sqrt(resvar),
+               modVariance = modVariance,
+               df = c(p, rdf), cov.unscaled = covar,
+               cov.scaled = covar * resvar,
+               info = object$info, niter = object$iterations,
+               stopmess = message,
+               par = param)
     
   ans$diffs <- object$diffs
   ans$data <- object$data
-  ans$additionalstats = CakeAdditionalStats(object$data)
+  ans$additionalstats <- CakeAdditionalStats(object$data)
   ans$seed <- object$seed
 
   ans$start <- object$start
   ans$fixed <- object$fixed
-  ans$errmin <- object$errmin 
+  ans$errmin <- object$errmin
+  ans$penalties <- object$penalties
   if(distimes){ 
     ans$distimes <- object$distimes
     ans$extraDT50 <- object$extraDT50
diff --git a/CakeModel.R b/CakeModel.R
index 3a32779..116ef6b 100644
--- a/CakeModel.R
+++ b/CakeModel.R
@@ -4,8 +4,8 @@
 # Portions Johannes Ranke, 2010
 # Contact: mkin-devel@lists.berlios.de
 
-# Modifications developed by Tessella for Syngenta: Copyright (C) 2011-2016 Syngenta
-# Tessella Project Reference: 6245, 7247, 8361, 7414
+# Modifications developed by Tessella for Syngenta: Copyright (C) 2011-2020 Syngenta
+# Tessella Project Reference: 6245, 7247, 8361, 7414, 10091
 
 #    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
@@ -22,6 +22,8 @@
 
 CakeModel <- function(..., use_of_ff = "max")
 {
+  # While spec is modified throughout this function, topology is added as-is to the returned model
+  topology <- list(...)
   spec <- list(...)
   obs_vars <- names(spec)
 
@@ -33,7 +35,7 @@ CakeModel <- function(..., use_of_ff = "max")
   parms <- vector()
   diffs <- vector()
   map <- list()
-  nlt = list()
+  nlt <- list()
 
   # Establish list of differential equations
   for (varname in obs_vars)
@@ -74,7 +76,7 @@ CakeModel <- function(..., use_of_ff = "max")
           ff <- vector()
           decline_term <- paste(nonlinear_term, "*", new_boxes[[1]])
         } else { # otherwise no decline term needed here
-          decline_term = "0" 
+          decline_term <- "0" 
         }
       } else {
         k_term <- paste("k", new_boxes[[1]], sep="_")
@@ -102,9 +104,9 @@ CakeModel <- function(..., use_of_ff = "max")
 
     # Construct and add DFOP term and add DFOP parameters if needed
     if(spec[[varname]]$type == "DFOP") {
-      k1_term = paste("k1", varname, sep="_")
-      k2_term = paste("k2", varname, sep="_")
-      g_term = paste("g", varname, sep="_")
+      k1_term <- paste("k1", varname, sep="_")
+      k2_term <- paste("k2", varname, sep="_")
+      g_term <- paste("g", varname, sep="_")
       
       new_parms <- c(k1_term, k2_term, g_term)
       ff <- vector()
@@ -121,7 +123,8 @@ CakeModel <- function(..., use_of_ff = "max")
         first_nonlinear_term <- paste(k1_term, new_boxes[[1]], sep=" * ")
         second_nonlinear_term <- paste(k2_term, new_boxes[[2]], sep=" * ")
         
-        overall_term_for_formation <- paste("(", g_term, " * ", first_nonlinear_term, " + (1 - ", g_term, ") *", second_nonlinear_term, ")")
+        # CU-150: formation has already been accounted for on the formation into the two sub-compartments; it does not need to be factored in again on the way out
+        overall_term_for_formation <- paste("(", first_nonlinear_term, " + ", second_nonlinear_term, ")")
         spec[[varname]]$nlt <- overall_term_for_formation
         
         new_diffs[[1]] <- paste(new_diffs[[1]], "-", first_nonlinear_term)
@@ -177,23 +180,23 @@ CakeModel <- function(..., use_of_ff = "max")
             # There is only one output, so no need for any flow fractions. Just add the whole flow from the parent
             formation_term <- spec[[varname]]$nlt
           } else {
-            fraction_to_target = paste("f", varname, "to", target, sep="_")
-            fraction_not_to_target = paste("(1 - ", fraction_to_target, ")", sep="")
+            fraction_to_target <- paste("f", varname, "to", target, sep="_")
+            fraction_not_to_target <- paste("(1 - ", fraction_to_target, ")", sep="")
             
             if(is.null(fraction_left)) {
-              fraction_really_to_target = fraction_to_target
-              fraction_left = fraction_not_to_target
+              fraction_really_to_target <- fraction_to_target
+              fraction_left <- fraction_not_to_target
             } else {
               # If this is the last output and there is no sink, it gets what's left
               if ( target==tail(to,1) && !spec[[varname]]$sink ) {
-                 fraction_really_to_target = fraction_left
+                 fraction_really_to_target <- fraction_left
               } else {
-                 fraction_really_to_target = fraction_to_target
-                 fraction_left = paste(fraction_left, " - ", fraction_to_target, sep="")
+                 fraction_really_to_target <- fraction_to_target
+                 fraction_left <- paste(fraction_left, " - ", fraction_to_target, sep="")
               }
             }
             
-            ff[target] = fraction_really_to_target
+            ff[target] <- fraction_really_to_target
             formation_term <- paste(ff[target], "*", spec[[varname]]$nlt)
             
             # Add the flow fraction parameter (if it exists)
@@ -218,9 +221,12 @@ CakeModel <- function(..., use_of_ff = "max")
     }
   }
   
-  model <- list(diffs = diffs, parms = parms, map = map, use_of_ff = use_of_ff)
+  model <- list(diffs = diffs, parms = parms, map = map, use_of_ff = use_of_ff, topology = topology)
 
-  if (exists("ff")) model$ff = ff
+  if (exists("ff")) {
+    model$ff <- ff
+  }
+  
   class(model) <- "mkinmod"
   invisible(model)
 }
diff --git a/CakeNullPlot.R b/CakeNullPlot.R
index e0823f2..52245ce 100644
--- a/CakeNullPlot.R
+++ b/CakeNullPlot.R
@@ -2,8 +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-2016 Syngenta
-# Tessella Project Reference: 6245, 7247, 8361, 7414
+# Modifications developed by Tessella for Syngenta: Copyright (C) 2011-2020 Syngenta
+# Tessella Project Reference: 6245, 7247, 8361, 7414, 10091
 
 #    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
@@ -19,7 +19,7 @@
 #    along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
 CakeNullPlot <- function(model, parms.ini, state.ini, observed, xlim = range(observed$time), 
-                         digits = max(3, getOption("digits") - 3), sannMaxIter = 10000, atol = 1e-6, ...)
+                         digits = max(3, getOption("digits") - 3), dfopDtMaxIter = 10000, atol = 1e-6, ...)
 {
   # Print the fixed values
   fixed <- data.frame(value = c(state.ini, parms.ini))
@@ -37,7 +37,7 @@ CakeNullPlot <- function(model, parms.ini, state.ini, observed, xlim = range(obs
   for (obs_var in obs_vars) {
     type = names(model$map[[obs_var]])[1]
     if(exists("DT50")) distimes[obs_var, ] = c(DT50, DT90)
-    distimes[obs_var, ] = CakeDT(type,obs_var,parms.all,sannMaxIter)
+    distimes[obs_var, ] = CakeDT(type,obs_var,parms.all,dfopDtMaxIter)
     extraDT50[obs_var, ] = CakeExtraDT(type, obs_var, parms.all)
   }
   
diff --git a/CakeOdeSolve.R b/CakeOdeSolve.R
new file mode 100644
index 0000000..263e57c
--- /dev/null
+++ b/CakeOdeSolve.R
@@ -0,0 +1,105 @@
+## -----------------------------------------------------------------------------
+## Ordinary Differential Equations Solver function.
+## -----------------------------------------------------------------------------
+# Some of the CAKE R modules are based on mkin.
+#
+# Modifications developed by Tessella for Syngenta, Copyright (C) 2011-2020 Syngenta
+# Tessella Project Reference: 6245, 7247, 8361, 7414, 10091
+#
+# The CAKE R modules are free software: you can
+# redistribute them and/or modify them under the
+# terms of the GNU General Public License as published by the Free Software
+# Foundation, either version 3 of the License, or (at your option) any later
+# version.
+#
+# This program is distributed in the hope that it will be useful, but WITHOUT
+# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+# FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
+# details.
+#
+# You should have received a copy of the GNU General Public License along with
+# this program. If not, see <http://www.gnu.org/licenses/> 
+
+
+# fit: Fit with initial (state) values and parameters for the ODE system.
+# outtimes: Time sequence for output.
+# solution: Whether to use analytical, eigenvectors, or general ode solver to solve the ODE system.
+# atol: The tolerance to apply to the ODE solver.
+CakeOdeSolve <- function(fit, outtimes, solution, atol)
+{
+  fixed <- fit$fixed$value
+  names(fixed) <- rownames(fit$fixed)
+  parms.all <- c(fit$par, fixed)
+  ininames <- c(
+    rownames(subset(fit$start, type == "state")),
+    rownames(subset(fit$fixed, type == "state")))
+  odeini <- parms.all[ininames]
+  names(odeini) <- gsub("_0$", "", names(odeini))
+  odenames <- c(
+    rownames(subset(fit$start, type == "deparm")),
+    rownames(subset(fit$fixed, type == "deparm")))
+  odeparms <- parms.all[odenames]
+  odeini <- AdjustOdeInitialValues(odeini, fit, odeparms)
+  
+  evalparse <- function(string)
+  {
+    eval(parse(text = string), as.list(c(odeparms, odeini)))
+  }
+  
+  odeResult <- numeric()
+  
+  if (solution == "analytical") 
+  {
+    parent.type = names(fit$map[[1]])[1]  
+    parent.name = names(fit$diffs)[[1]]
+    ode <- switch(parent.type,
+                  SFO = SFO.solution(outtimes, 
+                                     evalparse(parent.name),
+                                     evalparse(paste("k", parent.name, sep = "_"))),
+                  FOMC = FOMC.solution(outtimes,
+                                       evalparse(parent.name),
+                                       evalparse("alpha"), evalparse("beta")),
+                  DFOP = DFOP.solution(outtimes,
+                                       evalparse(parent.name),
+                                       evalparse(paste("k1", parent.name, sep = "_")), 
+                                       evalparse(paste("k2", parent.name, sep = "_")),
+                                       evalparse(paste("g",  parent.name, sep = "_"))),
+                  HS = HS.solution(outtimes,
+                                   evalparse(parent.name),
+                                   evalparse("k1"), evalparse("k2"),
+                                   evalparse("tb")),
+                  IORE = IORE.solution(outtimes,
+                                       evalparse(parent.name),
+                                       evalparse(paste("k", parent.name, sep = "_")),
+                                       evalparse("N"))
+         )
+    odeResult <- cbind(outtimes, ode)
+    dimnames(odeResult) <- list(outtimes, c("time", parent.name))
+  }
+  else if (solution == "eigen") 
+  {
+    coefmat.num <- matrix(sapply(as.vector(fit$coefmat), evalparse), 
+                          nrow = length(odeini))
+    e <- eigen(coefmat.num)
+    c <- solve(e$vectors, odeini)
+    f.out <- function(t) {
+      e$vectors %*% diag(exp(e$values * t), nrow = length(odeini)) %*% c
+    }
+    ode <- matrix(mapply(f.out, outtimes), 
+                  nrow = length(odeini), ncol = length(outtimes))
+    dimnames(ode) <- list(names(odeini), NULL)
+    odeResult <- cbind(time = outtimes, t(ode))
+  } 
+  else if (solution == "deSolve") 
+  {
+    odeResult <- ode(
+        y     = odeini,
+        times = outtimes,
+        func  = fit$mkindiff, 
+        parms = odeparms,
+        atol  = atol
+    )
+  }
+  
+  return(data.frame(odeResult))
+}
\ No newline at end of file
diff --git a/CakeOlsFit.R b/CakeOlsFit.R
index 2b8e736..e9b8f3c 100644
--- a/CakeOlsFit.R
+++ b/CakeOlsFit.R
@@ -7,18 +7,19 @@
 # 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-2016 Syngenta
-
+# Modifications developed by Tessella for Syngenta: Copyright (C) 2011-2020 Syngenta
+# Tessella Project Reference: 6245, 7247, 8361, 7414, 10091
+#
 #    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
 #    the Free Software Foundation, either version 3 of the License, or
 #    (at your option) any later version.
-# 
+#
 #    This program is distributed in the hope that it will be useful,
 #    but WITHOUT ANY WARRANTY; without even the implied warranty of
 #    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 #    GNU General Public License for more details.
-# 
+#
 #    You should have received a copy of the GNU General Public License
 #    along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
@@ -35,37 +36,37 @@
 # fixed_initials: A vector of compartments with fixed initial concentrations.
 # quiet: Whether the internal cost functions should execute more quietly than normal (less output).
 # atol: The tolerance to apply to the ODE solver.
-# sannMaxIter: The maximum number of iterations to apply to SANN processes.
-CakeOlsFit <- function(cake.model, 
+# 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)), 
-                       lower = 0, 
+                       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,
-                       sannMaxIter = 10000,
+                       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))
-  
+  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]
@@ -73,16 +74,16 @@ CakeOlsFit <- function(cake.model,
   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"
+    solution <- "analytical"
   } else {
-    solution = "deSolve"
+    solution <- "deSolve"
   }
-  
+
   # Create a function calculating the differentials specified by the model
   # if necessary
   if(solution == "deSolve") {
@@ -91,67 +92,71 @@ CakeOlsFit <- function(cake.model,
       diffs <- vector()
       for (box in mod_vars)
       {
-        diffname <- paste("d", box, sep="_")      
+        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
-  
-  out_predicted <- costFunctions$get.predicted()
-  
-  predicted_long <- wide_to_long(out_predicted, time = "time")
-  fit$predicted <- out_predicted
-  
+
   # 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$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$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
-  parms.all = c(fit$par, parms.fixed)
-  fit$distimes <- data.frame(DT50 = rep(NA, length(obs_vars)), DT90 = rep(NA, length(obs_vars)), 
+  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))
-  
+
   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,sannMaxIter)
-    fit$extraDT50[compartment.name, ] = CakeExtraDT(type, compartment.name, parms.all)
+    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_predicted, observed)
-  
+
+  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<-observed
   data$err<-rep(NA,length(data$time))
@@ -161,7 +166,8 @@ CakeOlsFit <- function(cake.model,
   data$variable <- ordered(data$variable, levels = obs_vars)
   fit$data <- data[order(data$variable, data$time), ]
   fit$atol <- atol
-  
-  class(fit) <- c("CakeFit", "mkinfit", "modFit") 
+  fit$topology <- cake.model$topology
+
+  class(fit) <- c("CakeFit", "mkinfit", "modFit")
   return(fit)
-}
\ No newline at end of file
+}
diff --git a/CakePenalties.R b/CakePenalties.R
index 5506b7f..b4f24d4 100644
--- a/CakePenalties.R
+++ b/CakePenalties.R
@@ -1,6 +1,7 @@
 # $Id$
 # Some of the CAKE R modules are based on mkin
-# CAKE (6245, 7247, 8361, 7414), by Tessella, for Syngenta: Copyright (C) 2011-2016 Syngenta
+# Modifications developed by Tessella for Syngenta: Copyright (C) 2011-2020 Syngenta
+# Tessella Project Reference: 6245, 7247, 8361, 7414, 10091
 #
 #    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
diff --git a/CakePlot.R b/CakePlot.R
index 1b80a4b..da69ea8 100755..100644
--- a/CakePlot.R
+++ b/CakePlot.R
@@ -1,8 +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-2016 Syngenta
-# Tessella Project Reference: 6245, 7247, 8361, 7414
+# Modifications developed by Tessella for Syngenta: Copyright (C) 2011-2020 Syngenta
+# Tessella Project Reference: 6245, 7247, 8361, 7414, 10091
 #
 #    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
@@ -19,21 +19,22 @@
 
 CakeDoPlot <- function(fit, xlim = range(fit$data$time), ...)
 {
+    original_fit <- fit
     t.map.names <- names(fit$map)
     metabolites <- grep("[A-Z]\\d", t.map.names, value = TRUE)
     t.map.rest  <- setdiff(t.map.names, metabolites)
     
     # Generate the normal graphs.
     final <- CakeSinglePlot(fit, xlim)
-    final_scaled <- final
     
-    if(length(metabolites) > 0){
-        for(i in 1:length(metabolites))
-        {
+    if (length(metabolites) > 0) {
+        for (i in 1:length(metabolites)) {
             metabolite <- metabolites[i]
             
-            if (names(fit$map[[metabolite]])[1] != "SFO"){
-                # We do not support these curves for models other than SFO (for now; see CU-79).
+            fit_method = names(fit$map[[metabolite]])[1]
+            if (fit_method != "SFO" && fit_method != "DFOP") {
+                # Only SFO and DFOP are current compatible with the process below
+                # so skip this metabolite.
                 next
             }
             
@@ -47,72 +48,49 @@ CakeDoPlot <- function(fit, xlim = range(fit$data$time), ...)
             
             metabolite0Parameter <- metabolite
             
-            if (!(metabolite %in% names(odeini))){
+            if (!(metabolite %in% names(odeini))) {
                 metabolite0Parameter <- paste0(metabolite, "_0")
             }
             
-            if (odeini[[metabolite0Parameter]] != 0){
-                # We do not support these curves where the initial concentration is non-zero (for now; see CU-79).
-                next
+            # If the metabolite is DFOP, we should find a k2_ decay variable,
+            # otherwise we can fall back to the SFO k_ variant
+            decay_var <- paste("k2", metabolite, sep="_")
+            if (!(decay_var %in% names(fit$par)) &&
+                !(decay_var %in% rownames(fit$fixed))) {
+                decay_var <- paste("k", metabolite, sep="_")
             }
             
-            decay_var <- paste("k", metabolite, sep="_")
-            
-            # calculate the new ffm (formation factor) and generate the two ffm scale charts
-            regex <- paste("f_(.+)_to", metabolite, sep="_")
-            decays = grep(regex, names(fit$par), value = TRUE)
-            
-            if (length(decays) != 1){
-                # We do not support these curves where there is formation from more than 1 compartment (for now; see CU-79).
-                next
-            }
-            
-            ffm_fitted <- sum(fit$par[decays])
-            normal <- final
-            ffm_scale <- NULL
-            
-            numeric_DT50 <- as.numeric(fit$distimes[metabolite,][["DT50"]])
-            
-            if (is.na(numeric_DT50)){
-                # We can't get anywhere without a numeric DT50.
-                next
-            }
-                
-            # Generate the curve for DT50=1000d and ffm as fitted.
-            if (decay_var %in% names(fit$par)){
-                k_new <- fit$par[decay_var]*numeric_DT50/1000;
-                fit$par[decay_var]<- k_new
-            }
-            else{
+            # Generate the curve for DT50=1000 and ffm as fitted.
+            if (decay_var %in% names(fit$par)) {
+                fit$par[decay_var] <- log(2)/1000
+            } else {
                 # If decay_var was fixed, need to modify it there.
-                k_new <- fit$fixed[decay_var,]["value"]*numeric_DT50/1000;
-                fit$fixed[decay_var,]["value"] <- k_new[decay_var,]
+                fit$fixed[decay_var,"value"] <- log(2)/1000
             }
             
-            dt50_1000_ffm_fitted <- CakeSinglePlot(fit, xlim)[metabolite]
-            
-            naming <- c(names(final), paste(metabolite, "DT50_1000_FFM_FITTED", sep="_"))
-            normal <- c(final, dt50_1000_ffm_fitted)
-            names(normal) <- naming
-            final <- normal
-            
-            # Generate the scaled FFM
-            if(ffm_fitted != 0)
-            {
-                ffm_scale <- 1 / ffm_fitted
-                final_scaled <- final[c("time", metabolite, paste(metabolite, "DT50_1000_FFM_FITTED", sep="_"))]
-                final_scaled[t.map.rest] <- NULL;
-                final_frame <- as.data.frame(final_scaled)
-                new_names <- c(paste(metabolite, "DT50_FITTED_FFM_1", sep="_"), paste(metabolite, "DT50_1000_FFM_1", sep="_"))
-                names(final_frame) <- c("time", new_names)
-                final_frame[new_names]<-final_frame[new_names]*ffm_scale;
-                
-                cat("<PLOT MODEL START>\n")
+            new_col_name <- paste(metabolite, "DT50_1000_FFM_FITTED", sep="_")
+            final[, new_col_name] <- CakeSinglePlot(fit, xlim)[metabolite]
             
-                write.table(final_frame, quote=FALSE)
+            ffm1_fit <- SetFormationFractionsToOne(fit, metabolite)
+            if (!is.null(ffm1_fit)) {
+                # It was possible to set the formation fractions to 1, so generate the curves
+                new_col_name <- paste(metabolite, "DT50_1000_FFM_1", sep="_")
+                final[, new_col_name] <- CakeSinglePlot(ffm1_fit, xlim)[metabolite]
                 
-                cat("<PLOT MODEL END>\n")
+                # Reset the k values back to their originals
+                # Then we can generate the curve for DT50 as fitted and ffm=1
+                if (decay_var %in% names(ffm1_fit$par)) {
+                    ffm1_fit$par[decay_var] <- original_fit$par[decay_var]
+                } else {
+                    ffm1_fit$fixed[decay_var,"value"] <- original_fit$fixed[decay_var,"value"]
+                }
+    
+                new_col_name <- paste(metabolite, "DT50_FITTED_FFM_1", sep="_")
+                final[, new_col_name] <- CakeSinglePlot(ffm1_fit, xlim)[metabolite]
             }
+            
+            # Finally reset the fit for the next metabolite
+            fit <- original_fit
         }
     }
     
@@ -125,98 +103,124 @@ CakeDoPlot <- function(fit, xlim = range(fit$data$time), ...)
     # View(final)
 }
 
-CakeSinglePlot <- function(fit, xlim = range(fit$data$time), scale_x = 1.0, ...)
+# Get the immediate parents of a metabolite in a topology
+GetParents <- function(topology, metabolite)
 {
-   solution = fit$solution
-   if ( is.null(solution) ) {
-      solution <- "deSolve"
-   }
-   atol = fit$atol
-   if ( is.null(atol) ) {
-      atol = 1.0e-6
-   }
-   
-  fixed <- fit$fixed$value
-  names(fixed) <- rownames(fit$fixed)
-  parms.all <- c(fit$par, fixed)
-  ininames <- c(
-    rownames(subset(fit$start, type == "state")),
-    rownames(subset(fit$fixed, type == "state")))
-  odeini <- parms.all[ininames]
-  names(odeini) <- gsub("_0$", "", names(odeini))
-  odenames <- c(
-    rownames(subset(fit$start, type == "deparm")),
-    rownames(subset(fit$fixed, type == "deparm")))
-  odeparms <- parms.all[odenames]
-  odeini <- AdjustOdeInitialValues(odeini, fit, odeparms)
+  parents <- NULL
+  
+  for (topo_metab in names(topology)) {
+    if (metabolite %in% topology[[topo_metab]]$to) {
+      parents <- append(parents, topo_metab)
+    }
+  }
+  
+  return (parents)
+}
 
-  outtimes <- seq(0, xlim[2], length.out=CAKE.plots.resolution) * scale_x
+# Get all the ancestors of a metabolite in a topology
+GetAncestors <- function(topology, metabolite, ancestors = NULL)
+{
+  for (parent in GetParents(topology, metabolite)) {
+    if (!(parent %in% ancestors)) {
+      ancestors <- GetAncestors(topology, parent, append(ancestors, parent))
+    }
+  }
+  
+  return (ancestors)
+}
 
-  # Solve the system
-  evalparse <- function(string)
-  {
-    eval(parse(text=string), as.list(c(odeparms, odeini)))
+# Set the formation leading to the metabolite to 1 from each parent
+# Returns the updated fit or NULL if the metabolite can't be updated for FFM=1
+SetFormationFractionsToOne <- function(fit, metabolite)
+{
+  parents <- GetParents(fit$topology, metabolite)
+
+  for (parent in parents) {
+    if (any(GetAncestors(fit$topology, parent) %in% parents)) {
+      # Reject metabolites where any parent is also an ancestor of another parent
+      return (NULL)
+    }
+    
+    # Try to find the formation fraction in the fitted and fixed parameters
+    ffm_var <- paste("f",  parent, "to", metabolite, sep="_")
+    updated_fit <- SetFfmVariable(fit, ffm_var, 1)
+    if (!is.null(updated_fit)) {
+      # Formation fraction was found and set to 1
+      fit <- updated_fit
+      next
+    }
+    
+    # The formation fraction for this metabolite is implicitly 1 minus all the others from this parent
+    # We need to iterate other formation fractions from the parent and set them to zero
+    for (other_metabolite in fit$topology[[parent]]$to) {
+      if (other_metabolite == metabolite) {
+        # Skip the metabolite itself because we know there isn't a formation fraction for it
+        next
+      }
+      
+      other_ffm_var <- paste("f", parent, "to", other_metabolite, sep="_")
+      updated_fit <- SetFfmVariable(fit, other_ffm_var, 0)
+      
+      # The updated fit shouldn't be able to be null because otherwise the topology was invalid
+      # or didn't match the parameters in the fit.  If it does, we'll throw an error.
+      if (is.null(updated_fit)) {
+        stop(paste("The updated fit was NULL. The formation fraction '", other_ffm_var, "' doesn't exist in the fit.", sep=""))
+      }
+      
+      fit <- updated_fit
+    }
   }
-  if (solution == "analytical") {
-    parent.type = names(fit$map[[1]])[1]  
-    parent.name = names(fit$diffs)[[1]]
-    o <- switch(parent.type,
-      SFO = SFO.solution(outtimes, 
-          evalparse(parent.name),
-          evalparse(paste("k", parent.name, sep="_"))),
-      FOMC = FOMC.solution(outtimes,
-          evalparse(parent.name),
-          evalparse("alpha"), evalparse("beta")),
-      DFOP = DFOP.solution(outtimes,
-          evalparse(parent.name),
-          evalparse(paste("k1", parent.name, sep="_")), 
-          evalparse(paste("k2", parent.name, sep="_")),
-          evalparse(paste("g", parent.name, sep="_"))),
-      HS = HS.solution(outtimes,
-          evalparse(parent.name),
-          evalparse("k1"), evalparse("k2"),
-          evalparse("tb")),
-      IORE = IORE.solution(outtimes,
-          evalparse(parent.name),
-          evalparse(paste("k", parent.name, sep="_")),
-          evalparse("N"))
-    )
-    out <- cbind(outtimes, o)
-    dimnames(out) <- list(outtimes, c("time", parent.name))
+  
+  return (fit)
+}
+
+# Set the value of ffm_var in the fit, looking in both fitted and fixed parameters
+# Returns the updated fit or NULL if the ffm_var couldn't be found
+SetFfmVariable <- function(fit, ffm_var, value)
+{
+  if (ffm_var %in% names(fit$par)) {
+    # The ffm_var was found as a fitted parameter
+    fit$par[ffm_var] <- value
+    return (fit)
   }
-  if (solution == "eigen") {
-    coefmat.num <- matrix(sapply(as.vector(fit$coefmat), evalparse), 
-      nrow = length(odeini))
-    e <- eigen(coefmat.num)
-    c <- solve(e$vectors, odeini)
-    f.out <- function(t) {
-      e$vectors %*% diag(exp(e$values * t), nrow=length(odeini)) %*% c
-    }
-    o <- matrix(mapply(f.out, outtimes), 
-      nrow = length(odeini), ncol = length(outtimes))
-    dimnames(o) <- list(names(odeini), NULL)
-    out <- cbind(time = outtimes, t(o))
-  } 
-  if (solution == "deSolve") {
-    out <- ode(
-      y = odeini,
-      times = outtimes,
-      func = fit$mkindiff, 
-      parms = odeparms,
-      atol = atol
-    )
+  
+  if (ffm_var %in% rownames(fit$fixed)) {
+    # The ffm_var was found as a fixed parameter
+    fit$fixed[ffm_var,"value"] <- value
+    return (fit)
+  }
+  
+  # ffm_var was not found in the fit
+  return (NULL)
+}
+
+CakeSinglePlot <- function(fit, xlim = range(fit$data$time), scale_x = 1.0, ...)
+{
+  solution <- fit$solution
+  if ( is.null(solution) ) {
+    solution <- "deSolve"
   }
   
-  out_transformed <- PostProcessOdeOutput(out, fit, atol)
+  atol <- fit$atol
+  if ( is.null(atol) ) {
+    atol <- 1.0e-5
+  }
+
+  outtimes <- seq(0, xlim[2], length.out=CAKE.plots.resolution) * scale_x
+
+  # Solve the system
+  odeResult <- CakeOdeSolve(fit, outtimes, solution, atol)
+  
+  odeResult_transformed <- PostProcessOdeOutput(odeResult, fit, atol)
   
   # Replace values that are incalculably small with 0.
   for (compartment.name in names(fit$map)) {
-    if (length(out_transformed[, compartment.name][!is.nan(out_transformed[, compartment.name])]) > 0) {
-      out_transformed[, compartment.name][is.nan(out_transformed[, compartment.name])] <- 0
+    if (length(odeResult_transformed[, compartment.name][!is.nan(odeResult_transformed[, compartment.name])]) > 0) {
+        odeResult_transformed[, compartment.name][is.nan(odeResult_transformed[, compartment.name])] <- 0
     }
     
-    out_transformed[, compartment.name][out_transformed[, compartment.name] < 0] <- 0
+    odeResult_transformed[, compartment.name][odeResult_transformed[, compartment.name] < 0] <- 0
   }
   
-  return(out_transformed)
+  return (odeResult_transformed)
 }
diff --git a/CakeSolutions.R b/CakeSolutions.R
index b78ab84..c9e5b88 100755..100644
--- a/CakeSolutions.R
+++ b/CakeSolutions.R
@@ -1,7 +1,7 @@
 # $Id$
 # Some of the CAKE R modules are based on mkin, 
-# Developed by Tessella Ltd for Syngenta: Copyright (C) 2011-2016 Syngenta
-# Tessella Project Reference: 6245, 7247, 8361, 7414
+# Developed by Tessella Ltd for Syngenta: Copyright (C) 2011-2020 Syngenta
+# Tessella Project Reference: 6245, 7247, 8361, 7414, 10091
 
 #    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
diff --git a/CakeSummary.R b/CakeSummary.R
index 70f01e6..5849a4a 100644
--- a/CakeSummary.R
+++ b/CakeSummary.R
@@ -3,8 +3,8 @@
 # and display the summary itself.
 
 # Parts modified from package mkin
-# Parts Tessella for Syngenta: Copyright (C) 2011-2016 Syngenta
-# Tessella Project Reference: 6245, 7247, 8361, 7414
+# Parts Tessella for Syngenta: Copyright (C) 2011-2020 Syngenta
+# Tessella Project Reference: 6245, 7247, 8361, 7414, 10091
 
 #    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
@@ -100,8 +100,8 @@ CakeFOMCBackCalculatedDT<-function(parms.all) {
 #  type - type of compartment
 #  obs_var - compartment name
 #  parms.all - list of parameters
-# sannMaxIter - the maximum amount of iterations to do for SANN
-CakeDT<-function(type, obs_var, parms.all, sannMaxIter) {    
+# dfopDtMaxIter - the maximum amount of iterations to do for DFOP DT calculation
+CakeDT<-function(type, obs_var, parms.all, dfopDtMaxIter) {    
     if (type == "SFO") {
       k_name = paste("k", obs_var, sep="_")
       k_tot = parms.all[k_name]
@@ -122,21 +122,19 @@ CakeDT<-function(type, obs_var, parms.all, sannMaxIter) {
         fraction <- g * exp( - k1 * t) + (1 - g) * exp( - k2 * t)
         (fraction - (1 - x/100))^2
       }
-
-      DTminBounds <- 0.001
       
-      # Determine a decent starting point for numerical iteration.  The latter two terms are also lower bounds for DT50.
-      DT50min <- max(0.001, (1 / k1) * log(g * 2), (1 / k2) * log((1 - g) * 2)) 
+      # Determine a decent starting point for numerical iteration.  These are lower bounds for DT50.
+      DT50min <- max((1 / k1) * log(g * 2), (1 / k2) * log((1 - g) * 2)) 
       
-      # Set the starting point for the numerical solver to a be a bit smaller than the computed minimum as the R optim function seems to get confused if the min is 
-      # greater than the answer it's trying to converge to (up to its level of accuracy).
+      # Set the starting point for the numerical solver to a be a bit smaller than the computed minimum.
       DT50Initial <- DT50min * 0.9
       
       # Results greater than 10,000 are not interesting.  The R optim routine also handles very large values unreliably and can claim to converge when it is nowhere near.
       if (DT50min > 10000){
           DT50 = ">10,000"
       }else{
-          DT50.temp <- optim(DT50Initial, f, method="SANN", control=list(fnscale=0.05, maxit=sannMaxIter), x = 50, lower = DTminBounds)
+          DT50.temp <- optim(DT50Initial, f, method="Nelder-Mead", control=list(maxit=dfopDtMaxIter, abstol=1E-10), x = 50)
+          
           DT50.o = DT50.temp$par
           DT50.converged = DT50.temp$convergence == 0
           DT50 = ifelse(!DT50.converged || DT50.o <= 0, NA, DT50.o)
@@ -146,17 +144,16 @@ CakeDT<-function(type, obs_var, parms.all, sannMaxIter) {
           }
       }
       
-      # Determine a decent starting point for numerical iteration.  The latter two terms are also lower bounds for DT90.
-      DT90min <- max(0.001, (1 / k1) * log(g * 10), (1 / k2) * log((1 - g) * 10))
+      # Determine a decent starting point for numerical iteration.  These are lower bounds for DT90.
+      DT90min <- max((1 / k1) * log(g * 10), (1 / k2) * log((1 - g) * 10))
       
-      # Set the starting point for the numerical solver to a be a bit smaller than the computed minimum as the R optim function seems to get confused if the min is 
-      # greater than the answer it's trying to converge to (up to its level of accuracy).
+      # Set the starting point for the numerical solver to a be a bit smaller than the computed minimum.
       DT90Initial <- DT90min * 0.9
       
       if (DT90min > 10000){
           DT90 = ">10,000"
       }else{
-          DT90.temp <- optim(DT90Initial, f, method="SANN", control=list(fnscale=0.05, maxit=sannMaxIter), x = 90, lower = DTminBounds)
+          DT90.temp <- optim(DT90Initial, f, method="Nelder-Mead", control=list(maxit=dfopDtMaxIter, abstol=1E-10), x = 90)
           DT90.o = DT90.temp$par
           DT90.converged = DT90.temp$convergence == 0
           DT90 = ifelse(!DT90.converged || DT90.o <= 0, NA, DT90.o)
@@ -321,7 +318,8 @@ summary.CakeFit <- function(object, data = TRUE, distimes = TRUE, halflives = TR
     dimnames(param) <- list(pnames, c("Estimate", "Std. Error",
                                     "t value", "Pr(>t)", "Lower CI (90%)", "Upper CI (90%)", "Lower CI (95%)", "Upper CI (95%)"))
                                     
-    ans <- list(residuals = object$residuals,
+    ans <- list(ssr = object$ssr,
+                residuals = object$residuals,
                 residualVariance = resvar,
                 sigma = sqrt(resvar),
                 modVariance = modVariance,
@@ -332,7 +330,8 @@ summary.CakeFit <- function(object, data = TRUE, distimes = TRUE, halflives = TR
                 par = param)
   } else {
     param <- cbind(param)  
-    ans <- list(residuals = object$residuals,
+    ans <- list(ssr = object$ssr,
+                residuals = object$residuals,
                 residualVariance = resvar,
                 sigma = sqrt(resvar),
                 modVariance = modVariance,
@@ -382,6 +381,8 @@ print.summary.CakeFit <- function(x, digits = max(3, getOption("digits") - 3), .
     cat("\nOptimised parameters:\n")
     printCoefmat(x$par, digits = digits, ...)
   }
+  
+  cat("\nSum of squared residuals:", format(signif(x$ssr, digits)))
 
   cat("\nResidual standard error:",
       format(signif(x$sigma, digits)), "on", rdf, "degrees of freedom\n")