#$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)
}
