Version 3.2v3.2

1 files changed, 84 insertions, 83 deletions

diff --git a/CakeNullPlot.R b/CakeNullPlot.R
index f5484b4..e0823f2 100644
--- a/CakeNullPlot.R
+++ b/CakeNullPlot.R
@@ -1,12 +1,11 @@
 #$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
+# Some of the CAKE R modules are based on mkin
 # Based on code in IRLSkinfit
-# Modifications developed by Tessella Plc for Syngenta: Copyright (C) 2011  Syngenta
-# Author: Rob Nelson (Tessella)
-# Tessella Project Reference: 6245
+# Modifications developed by Tessella for Syngenta: Copyright (C) 2011-2016 Syngenta
+# Tessella Project Reference: 6245, 7247, 8361, 7414
 
-#    This program is free software: you can redistribute it and/or modify
+#    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.
@@ -17,38 +16,57 @@
 #    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/>.”
+#    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, ...)
+CakeNullPlot <- function(model, parms.ini, state.ini, observed, xlim = range(observed$time), 
+                         digits = max(3, getOption("digits") - 3), sannMaxIter = 10000, atol = 1e-6, ...)
 {
   # Print the fixed values
   fixed <- data.frame(value = c(state.ini, parms.ini))
   fixed$type = c(rep("state", length(state.ini)), rep("deparm", length(parms.ini)))
   cat("\nFixed parameter values:\n")
   print(fixed)
+  parms.all = c(state.ini, parms.ini)
 
-  # Print disappearence times
+  # Print disappearance times
   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"))
+  extraDT50<- data.frame(k1 = rep(NA, length(obs_vars)), k2 = rep(NA, length(obs_vars)), row.names = obs_vars)
   
   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,sannMaxIter)
-   }
-
-   cat("\nEstimated disappearance times:\n")
-   print(distimes, digits=digits,...)
+    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)
+    extraDT50[obs_var, ] = CakeExtraDT(type, obs_var, parms.all)
+  }
+  
+  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
+  cat("\nHalf Lives for Extra k values:\n")
+  print(extraDT50, digits=digits,...)
+   
+  ioreRepDT = CakeIORERepresentativeDT("Parent", parms.all)
+  printIoreRepDT <- !is.null(ioreRepDT)
+  if (printIoreRepDT){
+    cat("\nRepresentative Half Life for IORE:", format(signif(ioreRepDT, digits)))
+  }
+   
+  cat("\n")
+  
+  fomcRepDT = CakeFOMCBackCalculatedDT(parms.all)
+  printFomcRepDT <- !is.null(fomcRepDT)
+  
+  if (printIoreRepDT){
+    cat("\nBack-calculated Half Life for FOMC:", format(signif(fomcRepDT, digits)))
+  }
+  
+  cat("\n")
+  
+  # Plot the model
 
-  outtimes <- seq(0, xlim[2], length.out=101)
+  outtimes <- seq(0, xlim[2], length.out=CAKE.plots.resolution)
 
   odeini <- state.ini
   odenames <- model$parms
@@ -60,75 +78,58 @@ CakeNullPlot <- function(model, parms.ini, state.ini, observed, xlim = range(obs
     eval(parse(text=string), as.list(c(odeparms, odeini)))
   }
   
-    mod_vars <- names(model$diffs)
-    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 = model$diffs[[box]])))
-        }
-        return(list(c(diffs)))
+  mod_vars <- names(model$diffs)
+  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 = model$diffs[[box]])))
     }
+    
+    return(list(c(diffs)))
+  }
+  
+  # Ensure initial state is at time 0
+  obstimes <- unique(c(0,observed$time))
+  
+  odeini <- AdjustOdeInitialValues(odeini, model, odeparms)
   
-   # Ensure initial state is at time 0
-   obstimes <- unique(c(0,observed$time))
-    out <- ode(
-      y = odeini,
-      times = obstimes,
-      func = mkindiff, 
-      parms = odeparms,
-    )
+  # Solve at observation times
+  out <- ode(y = odeini, times = obstimes, func = mkindiff, parms = odeparms, atol = atol)
     
-  # Output transformation for models with unobserved compartments like SFORB
-  out_transformed <- data.frame(time = out[,"time"])
-  for (var in names(model$map)) {
-    if(length(model$map[[var]]) == 1) {
-      out_transformed[var] <- out[, var]
-    } else {
-      out_transformed[var] <- rowSums(out[, model$map[[var]]])
-    }
-  }    
+  out_transformed <- PostProcessOdeOutput(out, model, atol)
+  
+  predicted_long <- wide_to_long(out_transformed,time="time")
   
-   predicted_long <- mkin_wide_to_long(out_transformed,time="time")
   # Calculate chi2 error levels according to FOCUS (2006)
-   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,...)
+  errmin<-CakeChi2(model, observed, predicted_long, obs_vars, c(), c(), state.ini, parms.ini, fixed)
+  cat("\nChi2 error levels in percent:\n")
+  errmin$err.min <- 100 * errmin$err.min
+  print(errmin, digits=digits,...)
 
+  # Fit to data
+  data<-observed
+  data$err<-rep(NA,length(data$time))
+  data<-merge(data, predicted_long, by=c("time","name"))
+  names(data)<-c("time", "variable", "observed","err-var", "predicted")
+  data$residual<-data$observed-data$predicted
+  data<-data[order(data$variable,data$time),]
+  cat("\nData:\n")
+  print(format(data, digits = digits, scientific = FALSE,...), row.names = FALSE)
   
-   # Fit to data
-   data<-observed
-   data$err<-rep(NA,length(data$time))
-   data<-merge(data, predicted_long, by=c("time","name"))
-   names(data)<-c("time", "variable", "observed","err-var", "predicted")
-   data$residual<-data$observed-data$predicted
-   data<-data[order(data$variable,data$time),]
-   cat("\nData:\n")
-   print(format(data, digits = digits, scientific = FALSE,...), row.names = FALSE)
-   
-  
-    out <- ode(
-      y = odeini,
-      times = outtimes,
-      func = mkindiff, 
-      parms = odeparms,
-    )
+  # Solve at output times
+  out <- ode(
+    y = odeini,
+    times = outtimes,
+    func = mkindiff, 
+    parms = odeparms)
     
-  # Output transformation for models with unobserved compartments like SFORB
-  out_transformed <- data.frame(time = out[,"time"])
-  for (var in names(model$map)) {
-    if(length(model$map[[var]]) == 1) {
-      out_transformed[var] <- out[, var]
-    } else {
-      out_transformed[var] <- rowSums(out[, model$map[[var]]])
-    }
-  }    
+  out_transformed <- PostProcessOdeOutput(out, model, atol)
    
-   cat("\n<PLOT MODEL START>\n")
-   print(out_transformed)
-   cat("<PLOT MODEL END>\n")
+  cat("\n<PLOT MODEL START>\n")
+  print(out_transformed)
+  cat("<PLOT MODEL END>\n")
 }