summaryrefslogtreecommitdiff
path: root/CakePlot.R
diff options
context:
space:
mode:
authorJohannes Ranke <jranke@uni-bremen.de>2017-10-18 11:28:39 +0200
committerJohannes Ranke <jranke@uni-bremen.de>2017-10-18 11:28:39 +0200
commit5b3daf393831acc4099e1bde3fe4527993529d74 (patch)
treea742cb6df0498fced89a7020467b99ad98fda468 /CakePlot.R
parent3d6b4b4b8293a4a4ab6f06805e1380600373796c (diff)
Version 3.2v3.2
Diffstat (limited to 'CakePlot.R')
-rwxr-xr-xCakePlot.R222
1 files changed, 222 insertions, 0 deletions
diff --git a/CakePlot.R b/CakePlot.R
new file mode 100755
index 0000000..1b80a4b
--- /dev/null
+++ b/CakePlot.R
@@ -0,0 +1,222 @@
+#$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
+#
+# 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/>.
+
+CakeDoPlot <- function(fit, xlim = range(fit$data$time), ...)
+{
+ 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))
+ {
+ 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).
+ next
+ }
+
+ 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]
+
+ metabolite0Parameter <- metabolite
+
+ 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
+ }
+
+ 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{
+ # 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,]
+ }
+
+ 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")
+
+ 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)
+}
+
+CakeSinglePlot <- function(fit, xlim = range(fit$data$time), scale_x = 1.0, ...)
+{
+ 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)
+
+ outtimes <- seq(0, xlim[2], length.out=CAKE.plots.resolution) * scale_x
+
+ # Solve the system
+ evalparse <- function(string)
+ {
+ eval(parse(text=string), as.list(c(odeparms, odeini)))
+ }
+ 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))
+ }
+ 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
+ )
+ }
+
+ out_transformed <- PostProcessOdeOutput(out, 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
+ }
+
+ out_transformed[, compartment.name][out_transformed[, compartment.name] < 0] <- 0
+ }
+
+ return(out_transformed)
+}

Contact - Imprint