Version 2.0v2.0

11 files changed, 548 insertions, 129 deletions

diff --git a/CakeCost.R b/CakeCost.R
index 1cc3981..40635f0 100644
--- a/CakeCost.R
+++ b/CakeCost.R
@@ -23,7 +23,7 @@
 # You should have received a copy of the GNU General Public License along with
 # this program. If not, see <http://www.gnu.org/licenses/># 
 #
-# $Id: CakeCost.R 216 2011-07-05 14:35:03Z nelr $
+# $Id$
 
 CakeCost <- function (model, obs, x = "time", y = NULL, err = NULL,
                      weight = "none", scaleVar = FALSE, cost = NULL, ...) {
@@ -180,6 +180,24 @@ CakeCost <- function (model, obs, x = "time", y = NULL, err = NULL,
     Res <- (ModVar- obsdat)
     res <- Res/Err
     resScaled <- res*Scale
+
+    # print("===========================");
+    # print("#Names[i]");
+    # print(Names[i]);
+    # print("xDat");
+    # print(xDat);
+    # print("obsdat");
+    # print(obsdat);
+    # print("ModVar");
+    # print(ModVar);
+    # print("1/Err");
+    # print(1/Err);
+    # print("Res");
+    # print(Res);
+    # print("res");
+    # print(res);
+    # print("===========================");
+    
     Residual <- rbind(Residual,
                       data.frame(
                         name   = Names[i],
@@ -191,13 +209,17 @@ CakeCost <- function (model, obs, x = "time", y = NULL, err = NULL,
                         res    = res))
 
     CostVar <- rbind(CostVar,
-                     data.frame(
-                       name           = Names[i],
-                       scale          = Scale,
-                       N              = length(Res),
-                       SSR.unweighted = sum(Res^2),
-                       SSR.unscaled   = sum(res^2),
-                       SSR            = sum(resScaled^2)))
+                  data.frame(
+                    name           = Names[i],
+                    scale          = Scale,
+                    N              = length(Res),
+                    SSR.unweighted = sum(Res^2),
+                    SSR.unscaled   = sum(res^2),
+                    SSR            = sum(resScaled^2)))
+                    
+    # print("************************")
+    # print(CostVar)
+    # print("************************")
   }  # end loop over all observed variables
 
   ## SSR
diff --git a/CakeInit.R b/CakeInit.R
index 3e609d5..0adab40 100644
--- a/CakeInit.R
+++ b/CakeInit.R
@@ -1,4 +1,4 @@
-#$Id: CakeInit.R 221 2011-11-01 14:12:45Z nelr $
+#$Id$
 
 # Purpose: Load the modules required by CAKE
 # Call with chdir=TRUE so it can load our source from the current directory
@@ -20,6 +20,7 @@
 #    along with this program.  If not, see <http://www.gnu.org/licenses/>.”
 
 library(mkin)
+library(Rsolnp)
 source("CakeModel.R")
 source("CakeOlsFit.R")
 source("CakeIrlsFit.R")
@@ -36,4 +37,4 @@ options(error=recover)
 options(warn=-1)
 
 # When running from C#, this must match the C# CAKE version
-CAKE.version<-"1.4"
+CAKE.version<-"2.0"
diff --git a/CakeIrlsFit.R b/CakeIrlsFit.R
index bf20572..05eff0a 100644
--- a/CakeIrlsFit.R
+++ b/CakeIrlsFit.R
@@ -1,9 +1,9 @@
-#$Id: CakeIrlsFit.R 216 2011-07-05 14:35:03Z nelr $
+#$Id$
 #
 # The CAKE R modules are based on mkin
 # Modifications developed by Tessella Plc for Syngenta: Copyright (C) 2011 Syngenta
-# Authors: Rob Nelson, Richard Smith
-# Tessella Project Reference: 6245
+# Authors: Rob Nelson, Richard Smith, Tamar Christina
+# Tessella Project Reference: 6245, 7247
  
 #    This program is free software: you can redistribute it and/or modify
 #    it under the terms of the GNU General Public License as published by
@@ -22,7 +22,8 @@ CakeIrlsFit <- function (mkinmod, observed, parms.ini = rep(0.1, length(mkinmod$
     state.ini = c(100, rep(0, length(mkinmod$diffs) - 1)), lower = 0, 
     upper = Inf, fixed_parms = NULL, fixed_initials = names(mkinmod$diffs)[-1], 
     plot = FALSE, quiet = FALSE, err = NULL, weight = "none", 
-    scaleVar = FALSE, atol=1e-6, control=list(),...) 
+    scaleVar = FALSE, atol=1e-6, sannMaxIter = 10000, control=list(),
+    useSolnp = FALSE, method='L-BFGS-B',...) 
 {
 
 ### This is a modification based on the "mkinfit" function.
@@ -35,21 +36,25 @@ CakeIrlsFit <- function (mkinmod, observed, parms.ini = rep(0.1, length(mkinmod$
     observed <- subset(observed, name %in% names(mkinmod$map))
     ERR <- rep(1,nrow(observed))
     observed <- cbind(observed,err=ERR)
+    flag <- 0
+    
     obs_vars = unique(as.character(observed$name))
     if (is.null(names(parms.ini))) 
         names(parms.ini) <- mkinmod$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 = mkinmod$diffs[[box]])))
+                                    eval(parse(text = mkinmod$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]
@@ -57,6 +62,7 @@ CakeIrlsFit <- function (mkinmod, observed, parms.ini = rep(0.1, length(mkinmod$
     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 = "_")
@@ -67,8 +73,37 @@ CakeIrlsFit <- function (mkinmod, observed, parms.ini = rep(0.1, length(mkinmod$
 	
     ############### Iteratively Reweighted Least Squares#############
     ## Start with no weighting
-    fit <- modFit(costFunctions$cost, c(state.ini.optim, parms.optim), lower = lower, 
-                      upper = upper,control=control,...)
+    
+    ## a prefitting step since this is usually the most effective method
+    if(useSolnp)
+    {
+        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)
+        #optimmethod <- method0
+
+        flag <- 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))
+            flag <- 0
+            if(class(a) == "try-error")
+            {
+                fit <- modFit(costFunctions$cost, c(state.ini.optim, parms.optim), lower = lower, upper = upper, method=method,control=control)
+            }
+        }
+    }else{        
+        fit <- modFit(costFunctions$cost, c(state.ini.optim, parms.optim), lower = lower, 
+                          upper = upper,control=control,...)
+    }
+    
+    ## print(fit$hessian)
 	
     if(length(control)==0)
       {
@@ -81,27 +116,109 @@ CakeIrlsFit <- function (mkinmod, observed, parms.ini = rep(0.1, length(mkinmod$
             control <- list(irls.control=irls.control)
           }
       }
+      
     irls.control <- control$irls.control
     maxIter <- irls.control$maxIter
     tol <- irls.control$tol
     ####
+    
+    if(length(mkinmod$map)==1 || useSolnp){
+        ## there is only one parent just do one iteration:
+        maxIter <- 0
+        
+        if(flag==1)## fit from 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
+        }
+        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(mod_vars))
     while(diffsigma>tol & niter<=maxIter)
       {      
+        if(flag==1 && useSolnp)## fit from 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
+        }
+
+        
         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
-        fit <- modFit(costFunctions$cost, fit$par, lower = lower, upper = upper, control=control, ...)
+        
+        if(useSolnp)
+        {
+            flag <- 1
+            a <- try(fit <- solnp(fit$par,fun=fn,LB=lower,UB=upper,control=control),silent=TRUE)
+
+            if(class(a) == "try-error")
+            {
+                flag <- 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{
+            fit <- modFit(costFunctions$cost, fit$par, lower = lower, upper = upper, control=control, ...)
+        }
+        
+        ## print(fit$hessian)
+                
         niter <- niter+1       
        
         ### If not converged, reweight and fit                
       }
+      
+    if(flag==1 && useSolnp){
+        ## 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
+    }
 	  
 	###########################################
     fit$mkindiff <- mkindiff
@@ -125,15 +242,19 @@ CakeIrlsFit <- function (mkinmod, observed, parms.ini = rep(0.1, length(mkinmod$
     fit$fixed$type = c(rep("state", length(state.ini.fixed)), 
         rep("deparm", length(parms.fixed)))
 
-    fit$errmin <- CakeChi2(observed, predicted_long, obs_vars, parms.optim, state.ini.optim)
+    fit$errmin <- CakeChi2(mkinmod, observed, predicted_long, obs_vars, parms.optim, state.ini.optim, state.ini, parms.ini)
     parms.all = c(fit$par, parms.fixed)
 	fit$penalties <- CakePenaltiesLong(parms.all, out_predicted, 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(DT50 = rep(NA, 2), row.names = c("k1", "k2"))   
+        
     for (obs_var in obs_vars) {
         type = names(mkinmod$map[[obs_var]])[1]
-        fit$distimes[obs_var, ] = CakeDT(type,obs_var,parms.all)
+        fit$distimes[obs_var, ] = CakeDT(type,obs_var,parms.all,sannMaxIter)
     }
+
+    fit$extraDT50[ ,c("DT50")] = CakeExtraDT(type,parms.all)
     
     data <- merge(observed, predicted_long, by = c("time", "name"))
     
@@ -142,6 +263,11 @@ CakeIrlsFit <- function (mkinmod, observed, parms.ini = rep(0.1, length(mkinmod$
     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
+    
     class(fit) <- c("CakeFit", "mkinfit", "modFit") 
     return(fit)
 }
diff --git a/CakeIrlsPlot.R b/CakeIrlsPlot.R
index 2de8910..5a7992f 100644
--- a/CakeIrlsPlot.R
+++ b/CakeIrlsPlot.R
@@ -1,4 +1,4 @@
-#$Id: CakeIrlsPlot.R 216 2011-07-05 14:35:03Z nelr $
+#$Id$
 # Generates fitted curves so the GUI can plot them
 # Author: Rob Nelson (Tessella)
 # Tessella Project Reference: 6245
@@ -20,6 +20,6 @@
 
 CakeIrlsPlot <- function(fit, xlim = range(fit$data$time), ...)
 {
-   CakeOlsPlot(fit, xlim, ...)
+   CakePlotInit(fit, xlim, ...)
 }
 
diff --git a/CakeMcmcFit.R b/CakeMcmcFit.R
index 209f904..f9a340e 100644
--- a/CakeMcmcFit.R
+++ b/CakeMcmcFit.R
@@ -1,9 +1,9 @@
-# $Id: CakeMcmcFit.R 216 2011-07-05 14:35:03Z nelr $
+# $Id$
 # The CAKE R modules are based on mkin, 
 # Based on mcmckinfit as modified by Bayer
 # Modifications developed by Tessella Plc for Syngenta: Copyright (C) 2011  Syngenta
-# Author: Rob Nelson
-# Tessella Project Reference: 6245
+# Author: Rob Nelson, Tamar Christina
+# Tessella Project Reference: 6245, 7247
 
 #    This program is free software: you can redistribute it and/or modify
 #    it under the terms of the GNU General Public License as published by
@@ -25,7 +25,9 @@ CakeMcmcFit <- function (mkinmod, observed, parms.ini = rep(0.1, length(mkinmod$
     scaleVar = FALSE, commonsigma=FALSE,jump = NULL,prior = NULL,
     wvar0=0.1,niter = 1000, 
     outputlength = niter, burninlength = 0, updatecov = niter, 
-    ntrydr = 1, drscale = NULL, verbose = TRUE,fitstart=TRUE,seed=NULL,atol=1e-6,...) 
+    ntrydr = 1, drscale = NULL, verbose = TRUE,fitstart=TRUE,seed=NULL,atol=1e-6,
+    sannMaxIter = 10000, control=list(),
+    useSolnp = FALSE, method='L-BFGS-B',...) 
 {
 
     NAind <-which(is.na(observed$value))
@@ -54,6 +56,7 @@ CakeMcmcFit <- function (mkinmod, observed, parms.ini = rep(0.1, length(mkinmod$
     state.ini.fixed <- state.ini[fixed_initials]
     optim_initials <- setdiff(names(state.ini), fixed_initials)
     state.ini.optim <- state.ini[optim_initials]
+    flag <- 0
     state.ini.optim.boxnames <- names(state.ini.optim)
     if (length(state.ini.optim) > 0) {
         names(state.ini.optim) <- paste(names(state.ini.optim), 
@@ -72,7 +75,7 @@ CakeMcmcFit <- function (mkinmod, observed, parms.ini = rep(0.1, length(mkinmod$
 		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
@@ -85,11 +88,44 @@ CakeMcmcFit <- function (mkinmod, observed, parms.ini = rep(0.1, length(mkinmod$
       {
      ## ############# Get Initial Paramtervalues   #############
     ## Start with no weighting
-    fit <- modFit(costFunctions$cost, c(state.ini.optim, parms.optim), lower = lower, 
-                      upper = upper,...)
-	
+    
+    if(useSolnp)
+    {
+        pnames=names(c(state.ini.optim, parms.optim))
+        fn <- function(P){
+            names(P) <- pnames
+            FF<<-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)
+        flag <- 1
+        optimmethod <- 'solnp'
+        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')
+            a <- try(fit <- modFit(costFunctions$cost, c(state.ini.optim, parms.optim), lower = lower, upper = upper, method='Port',control=control))
+            flag <- 0
+            
+            if(class(a) == "try-error")
+            {
+                fit <- modFit(costFunctions$cost, c(state.ini.optim, parms.optim), lower = lower, upper = upper, method=method,control=control)
+            }
+            ## now using submethod already
+        }
+    }
+    else
+    {
+        fit <- modFit(costFunctions$cost, c(state.ini.optim, parms.optim), lower = lower, 
+                          upper = upper,...)
+	}
+    
     if(commonsigma==TRUE){
       #observed$err <- rep(1,nrow(observed))
+      if(flag==1 && useSolnp)## fit from solnp
+      {
+          ## run a small loop with 'Marq' or some other method
+          fit <- modFit(costFunctions$cost, fit$par, lower = lower, upper = upper, method='Port',control=control)
+      }
       cov0 <- summary(fit)$cov.scaled*2.4^2/length(fit$par)
 	  costFunctions$set.calls(0); costFunctions$reset.best.cost()
       res <- modMCMC(costWithStatus,fit$par,...,jump=cov0,lower=lower,upper=upper,prior=prior,var0=var0,wvar0=wvar0,niter=niter,outputlength = outputlength, burninlength = burninlength, updatecov = updatecov,ntrydr=ntrydr,drscale=drscale,verbose=verbose)
@@ -111,6 +147,21 @@ CakeMcmcFit <- function (mkinmod, observed, parms.ini = rep(0.1, length(mkinmod$
       ERR <- err[as.character(observed$name)]
       observed$err <-ERR
 	  costFunctions$set.error(ERR)
+      if(flag==1 && useSolnp)## fit from 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
+      }
       olderr <- rep(1,length(mod_vars))
       diffsigma <- sum((err-olderr)^2)
       ## At least do one iteration step to get a weighted LS
@@ -166,11 +217,13 @@ CakeMcmcFit <- function (mkinmod, observed, parms.ini = rep(0.1, length(mkinmod$
    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(DT50 = rep(NA, 2), row.names = c("k1", "k2"))     
     for (obs_var in obs_vars) {
         type = names(mkinmod$map[[obs_var]])[1]
-        fit$distimes[obs_var, ] = CakeDT(type,obs_var,parms.all)
+        fit$distimes[obs_var, ] = CakeDT(type,obs_var,parms.all,sannMaxIter)
     }
 
+    fit$extraDT50[ ,c("DT50")] = CakeExtraDT(type,parms.all)
     
    # Fits to data
    odeini <- state.ini
@@ -208,7 +261,7 @@ CakeMcmcFit <- function (mkinmod, observed, parms.ini = rep(0.1, length(mkinmod$
 
    predicted_long <- mkin_wide_to_long(out_transformed,time="time")
    obs_vars = unique(as.character(observed$name))
-   fit$errmin <- CakeChi2(observed, predicted_long, obs_vars, parms.optim, state.ini.optim)
+   fit$errmin <- CakeChi2(mkinmod, observed, predicted_long, obs_vars, parms.optim, state.ini.optim, state.ini, parms.ini)
 
    data<-observed
    data$err<-rep(NA,length(data$time))
@@ -230,7 +283,7 @@ CakeMcmcFit <- function (mkinmod, observed, parms.ini = rep(0.1, length(mkinmod$
 
 
 # Summarise a fit
-summary.CakeMcmcFit <- function(object, data = TRUE, distimes = TRUE, ff = TRUE, cov = FALSE,...) {
+summary.CakeMcmcFit <- function(object, data = TRUE, distimes = TRUE, halflives = TRUE, ff = TRUE, cov = FALSE,...) {
   param  <- object$par
   pnames <- names(param)
   p      <- length(param)
@@ -255,12 +308,30 @@ summary.CakeMcmcFit <- function(object, data = TRUE, distimes = TRUE, ff = TRUE,
     warning(message)
   } else message <- "ok"
 
-   # Parameter fits
+    # 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.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
+    
+    # calculate the 95% confidence interval
     alpha <- 0.05
-    lci <- param + qt(alpha/2,rdf)*se
-    uci <- param + qt(1-alpha/2,rdf)*se
-   param <- cbind(param, se, tval, pt(tval, rdf, lower.tail = FALSE), lci, uci)
-   dimnames(param) <- list(pnames, c("Estimate", "Std. Error", "t value", "Pr(>t)", "Lower CI", "Upper CI"))
+    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
@@ -285,6 +356,7 @@ summary.CakeMcmcFit <- function(object, data = TRUE, distimes = TRUE, ff = TRUE,
   ans$fixed <- object$fixed
   ans$errmin <- object$errmin 
   if(distimes) ans$distimes <- object$distimes
+  if(halflives) ans$halflives <- object$halflives
   if(ff) ans$ff <- object$ff
   if(length(object$SFORB) != 0) ans$SFORB <- object$SFORB
   class(ans) <- c("summary.CakeFit","summary.mkinfit", "summary.modFit") 
diff --git a/CakeModel.R b/CakeModel.R
index a93d7af..ce6bc01 100644
--- a/CakeModel.R
+++ b/CakeModel.R
@@ -1,5 +1,5 @@
 # Was Id: mkinmod.R 71 2010-09-12 01:13:36Z jranke 
-# $Id: CakeModel.R 216 2011-07-05 14:35:03Z nelr $
+# $Id$
 
 # The CAKE R modules are based on mkin
 # Portions Johannes Ranke, 2010
@@ -22,12 +22,13 @@
 #    You should have received a copy of the GNU General Public License
 #    along with this program.  If not, see <http://www.gnu.org/licenses/>.”
 
-CakeModel <- function(...)
+CakeModel <- function(..., use_of_ff = "max")
 {
   spec <- list(...)
   obs_vars <- names(spec)
 
-  # The returned model will be a list of character vectors, containing
+  if (!use_of_ff %in% c("min", "max"))
+    stop("The use of formation fractions 'use_of_ff' can only be 'min' or 'max'")
   # differential equations, parameter names and a mapping from model variables
   # to observed variables. If possible, a matrix representation of the 
   # differential equations is included
@@ -69,21 +70,39 @@ CakeModel <- function(...)
 
     # Start a new differential equation for each new box
     new_diffs <- paste("d_", new_boxes, " =", sep="")
-
+    
+    # Get the name of the box(es) we are working on for the decline term(s)
+    box_1 = map[[varname]][[1]] # This is the only box unless type is SFORB
+    
     # Turn on sink if not specified otherwise
     if(is.null(spec[[varname]]$sink)) spec[[varname]]$sink <- TRUE
     
     #@@@ADD SFO k HERE !!!!!!!!!!!!!
     # Construct and add SFO term and add SFO parameters if needed
     if(spec[[varname]]$type == "SFO") {
+      if (use_of_ff == "min") { # Minimum use of formation fractions
+        if(spec[[varname]]$sink) {
       # From p. 53 of the FOCUS kinetics report
+      k_term <- paste("k", new_boxes[[1]], "sink", sep="_")
+      nonlinear_term <- paste(k_term, new_boxes[[1]], sep=" * ")
+      spec[[varname]]$nlt<-nonlinear_term
+      new_diffs[[1]] <- paste(new_diffs[[1]], "-", nonlinear_term)
+      new_parms <- k_term
+      ff <- vector()
+      decline_term <- paste(nonlinear_term, "*", new_boxes[[1]])
+    } else { # otherwise no decline term needed here
+      decline_term = "0" 
+    }
+  } else {
       k_term <- paste("k", new_boxes[[1]], sep="_")
       nonlinear_term <- paste(k_term, new_boxes[[1]], sep=" * ")
       spec[[varname]]$nlt<-nonlinear_term
       new_diffs[[1]] <- paste(new_diffs[[1]], "-", nonlinear_term)
       new_parms <- k_term
       ff <- vector()
-    } 
+    decline_term <- paste(nonlinear_term, "*", new_boxes[[1]])
+  }
+ } 
 
     # Construct and add FOMC term and add FOMC parameters if needed
     if(spec[[varname]]$type == "FOMC") {
@@ -225,33 +244,64 @@ CakeModel <- function(...)
       }
     }
   }
-  model <- list(diffs = diffs, parms = parms, map = map)
+  model <- list(diffs = diffs, parms = parms, map = map, use_of_ff = use_of_ff)
 
-  # Create coefficient matrix if appropriate
+  # Create coefficient matrix if appropriate#{{{
   if (mat) {
     boxes <- names(diffs)
     n <- length(boxes)
     m <- matrix(nrow=n, ncol=n, dimnames=list(boxes, boxes))
-    for (from in boxes) {
-      for (to in boxes) {
-        if (from == to) {
-          #@@@@ OMIT NEXT LINE? !!!!!
-          k.candidate = paste("k", from, c(boxes, "sink"), sep="_")
-          k.candidate = sub("free.*bound", "free_bound", k.candidate)
-          k.candidate = sub("bound.*free", "bound_free", k.candidate)
-          k.effective = intersect(model$parms, k.candidate)
-          m[from,to] = ifelse(length(k.effective) > 0,
-              paste("-", k.effective, collapse = " "), "0")
-        } else {
-          k.candidate = paste("k", from, to, sep="_")
-          k.candidate = sub("free.*bound", "free_bound", k.candidate)
-          k.candidate = sub("bound.*free", "bound_free", k.candidate)
-          k.effective = intersect(model$parms, k.candidate)
-          m[to, from] = ifelse(length(k.effective) > 0,
-              k.effective, "0")
+
+    if (use_of_ff == "min") { # Minimum use of formation fractions
+      for (from in boxes) {
+        for (to in boxes) {
+          if (from == to) { # diagonal elements
+            k.candidate = paste("k", from, c(boxes, "sink"), sep="_")
+            k.candidate = sub("free.*bound", "free_bound", k.candidate)
+            k.candidate = sub("bound.*free", "bound_free", k.candidate)
+            k.effective = intersect(model$parms, k.candidate)
+            m[from,to] = ifelse(length(k.effective) > 0,
+                paste("-", k.effective, collapse = " "), "0")
+
+          } else {          # off-diagonal elements
+            k.candidate = paste("k", from, to, sep="_")
+	    if (sub("_free$", "", from) == sub("_bound$", "", to)) {
+              k.candidate = paste("k", sub("_free$", "_free_bound", from), sep="_")
+	    }
+	    if (sub("_bound$", "", from) == sub("_free$", "", to)) {
+              k.candidate = paste("k", sub("_bound$", "_bound_free", from), sep="_")
+	    }
+            k.effective = intersect(model$parms, k.candidate)
+            m[to, from] = ifelse(length(k.effective) > 0,
+                k.effective, "0")
+          }
+        }
+      } 
+    } else {  # Use formation fractions where possible
+      for (from in boxes) {
+        for (to in boxes) {
+          if (from == to) { # diagonal elements
+            k.candidate = paste("k", from, sep="_")
+            m[from,to] = ifelse(k.candidate %in% model$parms,
+                paste("-", k.candidate), "0")
+            if(grepl("_free", from)) { # add transfer to bound compartment for SFORB
+              m[from,to] = paste(m[from,to], "-", paste("k", from, "bound", sep="_"))
+            }
+            if(grepl("_bound", from)) { # add backtransfer to free compartment for SFORB
+              m[from,to] = paste("- k", from, "free", sep="_")
+            }
+            m[from,to] = m[from,to]
+          } else {          # off-diagonal elements
+            f.candidate = paste("f", from, "to", to, sep="_")
+            k.candidate = paste("k", from, to, sep="_")
+	    # SFORB with maximum use of formation fractions not implemented, see above
+            m[to, from] = ifelse(f.candidate %in% model$parms,
+              paste(f.candidate, " * k_", from, sep=""), 
+              ifelse(k.candidate %in% model$parms, k.candidate, "0"))
+          }
         }
       }
-    }
+    } 
     model$coefmat <- m
   }
 
diff --git a/CakeNullPlot.R b/CakeNullPlot.R
index 3eacd67..f5484b4 100644
--- a/CakeNullPlot.R
+++ b/CakeNullPlot.R
@@ -1,4 +1,4 @@
-#$Id: CakeNullPlot.R 216 2011-07-05 14:35:03Z nelr $
+#$Id$
 # Generates model curves so the GUI can plot them. Used when all params are fixed so there was no fit
 # The CAKE R modules are based on mkin
 # Based on code in IRLSkinfit
@@ -19,7 +19,7 @@
 #    You should have received a copy of the GNU General Public License
 #    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), ...)
+CakeNullPlot <- function(model, parms.ini, state.ini, observed, xlim = range(observed$time), digits = max(3, getOption("digits") - 3), sannMaxIter = 10000, ...)
 {
   # Print the fixed values
   fixed <- data.frame(value = c(state.ini, parms.ini))
@@ -31,13 +31,20 @@ CakeNullPlot <- function(model, parms.ini, state.ini, observed, xlim = range(obs
   obs_vars = unique(as.character(observed$name))
   distimes <- data.frame(DT50 = rep(NA, length(obs_vars)), 
         DT90 = rep(NA, length(obs_vars)), row.names = obs_vars)
+  extraDT50<- data.frame(DT50 = rep(NA, 2), row.names = c("k1", "k2"))
+  
   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.ini)
+      distimes[obs_var, ] = CakeDT(type,obs_var,parms.ini,sannMaxIter)
    }
+
    cat("\nEstimated disappearance times:\n")
    print(distimes, digits=digits,...)
+   
+   extraDT50[ ,c("DT50")] = CakeExtraDT(type,parms.ini)
+   cat("\nHalf Lives for Extra k values:\n")
+   print(extraDT50, digits=digits,...)
 
    # Plot the model
 
@@ -45,7 +52,7 @@ CakeNullPlot <- function(model, parms.ini, state.ini, observed, xlim = range(obs
 
   odeini <- state.ini
   odenames <- model$parms
-  odeparms <- params.ini
+  odeparms <- parms.ini
 
   # Solve the system
   evalparse <- function(string)
@@ -86,7 +93,7 @@ CakeNullPlot <- function(model, parms.ini, state.ini, observed, xlim = range(obs
   
    predicted_long <- mkin_wide_to_long(out_transformed,time="time")
   # Calculate chi2 error levels according to FOCUS (2006)
-   errmin<-CakeChi2(observed, predicted_long, obs_vars, c(), c())
+   errmin<-CakeChi2(model, observed, predicted_long, obs_vars, c("auto"), c(), state.ini, c("auto"))
    cat("\nChi2 error levels in percent:\n")
    errmin$err.min <- 100 * errmin$err.min
    print(errmin, digits=digits,...)
diff --git a/CakeOlsFit.R b/CakeOlsFit.R
index 75ac471..626a595 100644
--- a/CakeOlsFit.R
+++ b/CakeOlsFit.R
@@ -7,11 +7,11 @@
 # from the FME package, v 1.1 by Soetart and Petzoldt, which was in turn
 # inspired by summary.nls.lm
 
-#$Id: CakeOlsFit.R 216 2011-07-05 14:35:03Z nelr $
+#$Id$
 # This version has been modified to expect SFO parameterised as k and flow fractions
 # Modifications developed by Tessella Plc for Syngenta: Copyright (C) 2011  Syngenta
-# Authors: Rob Nelson, Richard Smith
-# Tessella Project Reference: 6245
+# Authors: Rob Nelson, Richard Smith, Tamar Christina
+# Tessella Project Reference: 6245, 7247
 
 #    This program is free software: you can redistribute it and/or modify
 #    it under the terms of the GNU General Public License as published by
@@ -36,6 +36,8 @@ CakeOlsFit <- function(mkinmod, observed,
   plot = FALSE, quiet = FALSE,
   err = NULL, weight = "none", scaleVar = FALSE,
   atol = 1e-6,
+  sannMaxIter = 10000,
+  useSolnp = FALSE, method='L-BFGS-B',
   ...)
 {
   mod_vars <- names(mkinmod$diffs)
@@ -294,20 +296,23 @@ CakeOlsFit <- function(mkinmod, observed,
   fit$fixed$type = c(rep("state", length(state.ini.fixed)), rep("deparm", length(parms.fixed)))
 
 #  # Calculate chi2 error levels according to FOCUS (2006)
-  fit$errmin <- CakeChi2(observed, predicted_long, obs_vars, parms.optim, state.ini.optim)
+  fit$errmin <- CakeChi2(mkinmod, observed, predicted_long, obs_vars, parms.optim, state.ini.optim, state.ini, parms.ini)
 
   # 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(DT50 = rep(NA, 2), row.names = c("k1", "k2"))   
+        
   fit$ff <- vector()
   fit$SFORB <- vector()
   for (obs_var in obs_vars) {
     type = names(mkinmod$map[[obs_var]])[1]  
 
-    fit$distimes[obs_var, ] = CakeDT(type,obs_var,parms.all)
+    fit$distimes[obs_var, ] = CakeDT(type,obs_var,parms.all,sannMaxIter)
   }
 
+  fit$extraDT50[ ,c("DT50")] = CakeExtraDT(type,parms.all)  
   fit$penalties <- CakePenaltiesLong(parms.all, out_predicted, observed)
   
   # Collect observed, predicted and residuals
diff --git a/CakeOlsPlot.R b/CakeOlsPlot.R
index 199cc28..e92742c 100644
--- a/CakeOlsPlot.R
+++ b/CakeOlsPlot.R
@@ -1,4 +1,4 @@
-#$Id: CakeOlsPlot.R 216 2011-07-05 14:35:03Z nelr $
+#$Id$
 # Generates fitted curves so the GUI can plot them
 # Based on code in IRLSkinfit
 # Author: Rob Nelson (Tessella)
@@ -18,10 +18,70 @@
 #    You should have received a copy of the GNU General Public License
 #    along with this program.  If not, see <http://www.gnu.org/licenses/>.”
 
-CakeOlsPlot <- function(fit, xlim = range(fit$data$time), ...)
+CakePlotInit <- function(fit, xlim = range(fit$data$time), ...)
 {
-   cat("<PLOT MODEL START>\n")
+    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 <- CakeOlsPlot(fit, xlim)
+    final_scaled <- final
+    
+    if(length(metabolites) > 0){
+        for(i in 1:length(metabolites))
+        {
+            metabolite <- metabolites[i]
+            decay_var  <- paste("k", metabolite, sep="_")
+            
+            # calculate the new ffm and generate the two ffm scale charts
+            regex <- paste("f_(.+)_to", metabolite, sep="_")
+            decays = grep(regex, names(fit$par), value = TRUE)
+            ffm_fitted <- sum(fit$par[decays])
+            normal <- final
+            ffm_scale <- NULL
+                
+            # Generate the DT50=1000d and ffm as fitted.
+            k_new <- fit$par[decay_var]*fit$distimes[metabolite,]["DT50"]/1000;
+            fit$par[decay_var]<- k_new[metabolite,]
+            dt50_1000_ffm_fitted <- CakeOlsPlot(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")
+            
+                write.table(final_frame, quote=FALSE)
+                
+                cat("<PLOT MODEL END>\n")
+            }
+        }
+    }
+    
+    cat("<PLOT MODEL START>\n")
+    
+    write.table(final, quote=FALSE)
+    
+    cat("<PLOT MODEL END>\n")
+
+    # View(final)
+}
 
+CakeOlsPlot <- function(fit, xlim = range(fit$data$time), scale_x = 1.0, ...)
+{
    solution = fit$solution
    if ( is.null(solution) ) {
       solution <- "deSolve"
@@ -40,7 +100,7 @@ CakeOlsPlot <- function(fit, xlim = range(fit$data$time), ...)
   odeini <- parms.all[ininames]
   names(odeini) <- names(fit$diffs)
 
-  outtimes <- seq(0, xlim[2], length.out=101)
+  outtimes <- seq(0, xlim[2], length.out=101) * scale_x
 
   odenames <- c(
     rownames(subset(fit$start, type == "deparm")),
@@ -111,7 +171,5 @@ CakeOlsPlot <- function(fit, xlim = range(fit$data$time), ...)
       out_transformed[var] <- rowSums(out[, fit$map[[var]]])
     }
   }
-  print(out_transformed)
-
-  cat("<PLOT MODEL END>\n")
+  return(out_transformed)
 }
diff --git a/CakePenalties.R b/CakePenalties.R
index 3b2a0df..1feffbe 100644
--- a/CakePenalties.R
+++ b/CakePenalties.R
@@ -1,4 +1,4 @@
-# $Id: CakePenalties.R 216 2011-07-05 14:35:03Z nelr $
+# $Id$
 # The CAKE R modules are based on mkin
 # CAKE (6245), by Tessella, for Syngenta: Copyright (C) 2011  Syngenta
 #
diff --git a/CakeSummary.R b/CakeSummary.R
index ec42e75..6c7c661 100644
--- a/CakeSummary.R
+++ b/CakeSummary.R
@@ -1,10 +1,10 @@
-# $Id: CakeSummary.R 216 2011-07-05 14:35:03Z nelr $
+# $Id$
 # CakeSummary: Functions to calculate statistics used in the summary,
 # and display the summary itself.
 
 # Parts modified from package mkin
-# Parts Tessella (Rob Nelson/Richard Smith) for Syngenta: Copyright (C) 2011  Syngenta
-# Tessella Project Reference: 6245
+# Parts Tessella (Rob Nelson/Richard Smith/Tamar Christina) for Syngenta: Copyright (C) 2011  Syngenta
+# Tessella Project Reference: 6245, 7247
 
 #    This program is free software: you can redistribute it and/or modify
 #    it under the terms of the GNU General Public License as published by
@@ -19,12 +19,32 @@
 #    You should have received a copy of the GNU General Public License
 #    along with this program.  If not, see <http://www.gnu.org/licenses/>.”
 
+# Calculate the extra decay times for the other k values in DFOP
+# and HS
+# Arguments
+#  type - type of compartment
+#  parms.all - list of parameters
+CakeExtraDT<-function(type, parms.all) {  
+    DT50k1 = NA    
+    DT50k2 = NA    
+    if (!is.na(parms.all["k1"]) && !is.na(parms.all["k2"])) {
+      k1 = parms.all["k1"]
+      k2 = parms.all["k2"]
+      DT50k1 = log(2)/k1
+      DT50k2 = log(2)/k2
+    }
+    ret<-c(DT50k1, DT50k2)
+    names(ret)<-c("k1", "k2")
+    ret
+}
+
 # Calculate the decay times for a compartment
 # Arguments
 #  type - type of compartment
 #  obs_var - compartment name
 #  parms.all - list of parameters
-CakeDT<-function(type, obs_var, parms.all) {
+# sannMaxIter - the maximum amount of iterations to do for SANN
+CakeDT<-function(type, obs_var, parms.all, sannMaxIter) {    
     if (type == "SFO") {
       # Changed to new parameterisation
       #k_names = grep(paste("k", obs_var, sep="_"), names(parms.all), value=TRUE)
@@ -47,11 +67,17 @@ CakeDT<-function(type, obs_var, parms.all) {
         fraction <- g * exp( - k1 * t) + (1 - g) * exp( - k2 * t)
         (fraction - (1 - x/100))^2
       }
-      DTmax <- 1000
-      DT50.o <- optimize(f, c(0.001, DTmax), x=50)$minimum
-      DT50 = ifelse(DTmax - DT50.o < 0.1, NA, DT50.o)
-      DT90.o <- optimize(f, c(0.001, DTmax), x=90)$minimum
-      DT90 = ifelse(DTmax - DT90.o < 0.1, NA, DT90.o)
+
+      DTminBounds <- 0.001
+      DTminInitial <- 1
+      DT50.temp <- optim(DTminInitial, f, method="SANN", control=list(fnscale=0.05, maxit=sannMaxIter), x = 50, lower = DTminBounds)
+      DT50.o = DT50.temp$par
+      DT50.converged = DT50.temp$convergence == 0
+      DT50 = ifelse(!DT50.converged, NA, DT50.o)
+      DT90.temp <- optim(DTminInitial, f, method="SANN", control=list(fnscale=0.05, maxit=sannMaxIter), x = 90, lower = DTminBounds)
+      DT90.o = DT90.temp$par
+      DT90.converged = DT90.temp$convergence == 0
+      DT90 = ifelse(!DT90.converged, NA, DT90.o)
     }
     if (type == "HS") {
       k1 = parms.all["k1"]
@@ -84,12 +110,16 @@ CakeDT<-function(type, obs_var, parms.all) {
         ((k_12 + k_21 - b2)/(b1 - b2)) * exp(-b2 * t)
       }
       f_50 <- function(t) (SFORB_fraction(t) - 0.5)^2
-      max_DT <- 1000
-      DT50.o <- optimize(f_50, c(0.01, max_DT))$minimum
-      if (abs(DT50.o - max_DT) < 0.01) DT50 = NA else DT50 = DT50.o
+      DTmin <- 0.001
+      DT50.temp <- optim(DTmin, f_50, method="SANN", control=list(fnscale=0.05, maxit=sannMaxIter), x = 50, lower = DTmin)
+      DT50.o = DT50.temp$par
+      DT50.converged = DT50.temp$convergence == 0
+      if (!DT50.converged) DT50 = NA else DT50 = DT50.o
       f_90 <- function(t) (SFORB_fraction(t) - 0.1)^2
-      DT90.o <- optimize(f_90, c(0.01, 1000))$minimum
-      if (abs(DT90.o - max_DT) < 0.01) DT90 = NA else DT90 = DT90.o
+      DT90.temp <- optim(DTmin, f_90, method="SANN", control=list(fnscale=0.05, maxit=sannMaxIter), x = 90, lower = DTmin)
+      DT90.o = DT90.temp$par
+      DT90.converged = DT90.temp$convergence == 0
+      if (!DT90.converged) DT90 = NA else DT90 = DT90.o
     }
     ret<-c(DT50, DT90)
     names(ret)<-c("DT50","DT90")
@@ -104,32 +134,48 @@ CakeDT<-function(type, obs_var, parms.all) {
 #  parms.optim - list of fitted parameters
 #  state,ini.optim - list of fitted initial values
 #
-CakeChi2<-function(observed, predicted_long, obs_vars, parms.optim, state.ini.optim) {
+CakeChi2<-function(mkinmod, observed, predicted_long, obs_vars, parms.optim, state.ini.optim, state.ini, parms.ini) {
   # Calculate chi2 error levels according to FOCUS (2006)
-  means <- aggregate(value ~ time + name, data = observed, mean, na.rm=TRUE)
-  errdata <- merge(means, predicted_long, by = c("time", "name"), suffixes = c("_mean", "_pred"))
-  errdata <- errdata[order(errdata$time, errdata$name), ]
-  errmin.overall <- mkinerrmin(errdata, length(parms.optim) + length(state.ini.optim))
-
-  errmin <- data.frame(err.min = errmin.overall$err.min, 
-    n.optim = errmin.overall$n.optim, df = errmin.overall$df)
-  rownames(errmin) <- "All data"
-  for (obs_var in obs_vars)
-  {
-    errdata.var <- subset(errdata, name == obs_var)
-    n.k.optim <- length(grep(paste("k", obs_var, sep="_"), names(parms.optim)))
-    n.initials.optim <- length(grep(paste(obs_var, ".*", "_0", sep=""), names(state.ini.optim)))
-    n.optim <- n.k.optim + n.initials.optim
-    if ("alpha" %in% names(parms.optim)) n.optim <- n.optim + 1
-    if ("beta" %in% names(parms.optim)) n.optim <- n.optim + 1
-    if ("k1" %in% names(parms.optim)) n.optim <- n.optim + 1
-    if ("k2" %in% names(parms.optim)) n.optim <- n.optim + 1
-    if ("g" %in% names(parms.optim)) n.optim <- n.optim + 1
-    if ("tb" %in% names(parms.optim)) n.optim <- n.optim + 1
-    errmin.tmp <- mkinerrmin(errdata.var, n.optim)
-    errmin[obs_var, c("err.min", "n.optim", "df")] <- errmin.tmp
+    
+  if(getRversion() >= '2.14.0'){
+    # You would usually call mkinfit to fit the data, however it seems that we've already fit the data
+    # bundle <- mkinfit(mkinmod, observed, parms.ini = parms.optim, state.ini = state.ini) #, err='err') 
+    # Somewhere else. Calling mkinfit may result in an error during fitting. So Instead we just create
+    # the values that the new version of mkinerrmin expects. There is no class check so it's fine.
+    bundle<-vector()
+    bundle$predicted <- predicted_long
+    bundle$observed <- observed
+    bundle$obs_vars <- obs_vars
+    bundle$par <- c(parms.optim, state.ini.optim)
+    errmin.overall <- mkinerrmin(bundle)
+    
+    errmin.overall
+  } else {
+    means <- aggregate(value ~ time + name, data = observed, mean, na.rm=TRUE)
+    errdata <- merge(means, predicted_long, by = c("time", "name"), suffixes = c("_mean", "_pred"))
+    errdata <- errdata[order(errdata$time, errdata$name), ]
+    errmin.overall <- mkinerrmin(errdata, length(parms.optim) + length(state.ini.optim))
+    
+    errmin <- data.frame(err.min = errmin.overall$err.min, 
+      n.optim = errmin.overall$n.optim, df = errmin.overall$df)
+    rownames(errmin) <- "All data"
+    for (obs_var in obs_vars)
+    {
+      errdata.var <- subset(errdata, name == obs_var)
+      n.k.optim <- length(grep(paste("k", obs_var, sep="_"), names(parms.optim)))
+      n.initials.optim <- length(grep(paste(obs_var, ".*", "_0", sep=""), names(state.ini.optim)))
+      n.optim <- n.k.optim + n.initials.optim
+      if ("alpha" %in% names(parms.optim)) n.optim <- n.optim + 1
+      if ("beta" %in% names(parms.optim)) n.optim <- n.optim + 1
+      if ("k1" %in% names(parms.optim)) n.optim <- n.optim + 1
+      if ("k2" %in% names(parms.optim)) n.optim <- n.optim + 1
+      if ("g" %in% names(parms.optim)) n.optim <- n.optim + 1
+      if ("tb" %in% names(parms.optim)) n.optim <- n.optim + 1
+      errmin.tmp <- mkinerrmin(errdata.var, n.optim)
+      errmin[obs_var, c("err.min", "n.optim", "df")] <- errmin.tmp
+    }
+    errmin
   }
-  errmin
 }
 
 # Compute additional statistics (rē/efficiency) of obs v pred for each compartment and the whole problem
@@ -152,15 +198,20 @@ CakeAdditionalStats<-function(obsvspred){
 }
 
 # Summarise a fit - used for CakeIrlsFit and CakeOlsFit
-summary.CakeFit <- function(object, data = TRUE, distimes = TRUE, ff = TRUE, cov = FALSE,...) {
+summary.CakeFit <- 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
+  # covar  <- try(ginv(0.5*object$hessian))   # unscaled covariance, calculate using generic inversion
 
   rdf    <- object$df.residual
   resvar <- object$ssr / rdf
   modVariance <- object$ssr / length(object$residuals)
+    
+  # if(!is.numeric(covar)){
+   # covar  <- try(ginv(0.5*object$hessian))
+  # }
   
   if (!is.numeric(covar)) {
 #    message <- "Cannot estimate covariance; system is singular"
@@ -181,13 +232,33 @@ summary.CakeFit <- function(object, data = TRUE, distimes = TRUE, ff = TRUE, cov
     warning(message)
   } else message <- "ok"
 
-  if(cov && !is.null(covar)) {
+  if(cov && !is.null(covar)) { 
+
+    # 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.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
+    
+    # calculate the 95% confidence interval
     alpha <- 0.05
-    lci <- param + qt(alpha/2,rdf)*se
-    uci <- param + qt(1-alpha/2,rdf)*se
-    param <- cbind(param, se, tval, pt(tval, rdf, lower.tail = FALSE), lci, uci)
+    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", "Upper CI"))
+                                    "t value", "Pr(>t)", "Lower CI (90%)", "Upper CI (90%)", "Lower CI (95%)", "Upper CI (95%)"))
+                                    
     ans <- list(residuals = object$residuals,
                 residualVariance = resvar,
                 sigma = sqrt(resvar),
@@ -219,14 +290,15 @@ summary.CakeFit <- function(object, data = TRUE, distimes = TRUE, ff = TRUE, cov
   ans$errmin <- object$errmin 
   ans$penalties <- object$penalties
   if(distimes) ans$distimes <- object$distimes
+  if(halflives) ans$extraDT50 <- object$extraDT50
   if(ff) ans$ff <- object$ff
   if(length(object$SFORB) != 0) ans$SFORB <- object$SFORB
-  class(ans) <- c("summary.CakeFit","summary.mkinfit", "summary.modFit") 
+  class(ans) <- c("summary.CakeFit","summary.mkinfit", "summary.cakeModFit") 
   return(ans)  
 }
 
 # Print a summary. Used for CakeIrlsFit, CakeOlsFit and CakeMcmcFit
-# Expanded from print.summary.modFit
+# Expanded from print.summary.cakeModFit
 print.summary.CakeFit <- function(x, digits = max(3, getOption("digits") - 3), ...) {
   cat("\nEquations:\n")
   print(noquote(as.character(x[["diffs"]])))
@@ -239,7 +311,7 @@ print.summary.CakeFit <- function(x, digits = max(3, getOption("digits") - 3), .
   cat("\nFixed parameter values:\n")
   if(length(x$fixed$value) == 0) cat("None\n")
   else print(x$fixed)
-  
+
   if ( !is.null(x$par) ) {
     cat("\nOptimised parameters:\n")
     printCoefmat(x$par, digits = digits, ...)
@@ -256,7 +328,13 @@ print.summary.CakeFit <- function(x, digits = max(3, getOption("digits") - 3), .
   if(printdistimes){
     cat("\nEstimated disappearance times:\n")
     print(x$distimes, digits=digits,...)
-  }    
+  }  
+  
+  printextraDT50 <- !is.null(x$extraDT50)
+  if(printextraDT50){
+    cat("\nHalf Lives for Extra k values:\n")
+    print(x$extraDT50, digits=digits,...)
+  }  
 
   printff <- !is.null(x$ff)
   if(printff){