diff options
Diffstat (limited to 'CakePlot.R')
| -rwxr-xr-x | CakePlot.R | 222 |
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) +} |