summaryrefslogtreecommitdiff
path: root/CakeNullPlot.R
diff options
context:
space:
mode:
Diffstat (limited to 'CakeNullPlot.R')
-rw-r--r--CakeNullPlot.R167
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")
}

Contact - Imprint