1 files changed, 127 insertions, 0 deletions
diff --git a/CakeNullPlot.R b/CakeNullPlot.R
new file mode 100644
index 0000000..3eacd67
--- /dev/null
+++ b/CakeNullPlot.R
@@ -0,0 +1,127 @@
+#$Id: CakeNullPlot.R 216 2011-07-05 14:35:03Z nelr $
+# 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
+# Modifications developed by Tessella Plc for Syngenta: Copyright (C) 2011  Syngenta
+# Author: Rob Nelson (Tessella)
+# Tessella Project Reference: 6245
+
+#    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
+#    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/>.”
+
+CakeNullPlot <- function(model, parms.ini, state.ini, observed, xlim = range(observed$time), digits = max(3, getOption("digits") - 3), ...)
+{
+  # 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)
+
+  # Print disappearence 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)
+  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)
+   }
+   cat("\nEstimated disappearance times:\n")
+   print(distimes, digits=digits,...)
+
+   # Plot the model
+
+  outtimes <- seq(0, xlim[2], length.out=101)
+
+  odeini <- state.ini
+  odenames <- model$parms
+  odeparms <- params.ini
+
+  # Solve the system
+  evalparse <- function(string)
+  {
+    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)))
+    }
+  
+   # Ensure initial state is at time 0
+   obstimes <- unique(c(0,observed$time))
+    out <- ode(
+      y = odeini,
+      times = obstimes,
+      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]]])
+    }
+  }    
+  
+   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())
+   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)
+   
+  
+    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]]])
+    }
+  }    
+   
+   cat("\n<PLOT MODEL START>\n")
+   print(out_transformed)
+   cat("<PLOT MODEL END>\n")
+}
+