#$Id$
# Generates fitted curves so the GUI can plot them
# Based on code in IRLSkinfit
# Modifications developed by Hybrid Intelligence (formerly Tessella), part of
# Capgemini Engineering, for Syngenta, Copyright (C) 2011-2022 Syngenta
# Tessella Project Reference: 6245, 7247, 8361, 7414, 10091
#
# 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), ...)
{
original_fit <- fit
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)
if (length(metabolites) > 0) {
for (i in 1:length(metabolites)) {
metabolite <- metabolites[i]
fit_method = names(fit$map[[metabolite]])[1]
if (fit_method != "SFO" && fit_method != "DFOP") {
# Only SFO and DFOP are current compatible with the process below
# so skip this metabolite.
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 the metabolite is DFOP, we should find a k2_ decay variable,
# otherwise we can fall back to the SFO k_ variant
decay_var <- paste("k2", metabolite, sep="_")
if (!(decay_var %in% names(fit$par)) &&
!(decay_var %in% rownames(fit$fixed))) {
decay_var <- paste("k", metabolite, sep="_")
}
# Generate the curve for DT50=1000 and ffm as fitted.
if (decay_var %in% names(fit$par)) {
fit$par[decay_var] <- log(2)/1000
} else {
# If decay_var was fixed, need to modify it there.
fit$fixed[decay_var,"value"] <- log(2)/1000
}
new_col_name <- paste(metabolite, "DT50_1000_FFM_FITTED", sep="_")
final[, new_col_name] <- CakeSinglePlot(fit, xlim)[metabolite]
ffm1_fit <- SetFormationFractionsToOne(fit, metabolite)
if (!is.null(ffm1_fit)) {
# It was possible to set the formation fractions to 1, so generate the curves
new_col_name <- paste(metabolite, "DT50_1000_FFM_1", sep="_")
final[, new_col_name] <- CakeSinglePlot(ffm1_fit, xlim)[metabolite]
# Reset the k values back to their originals
# Then we can generate the curve for DT50 as fitted and ffm=1
if (decay_var %in% names(ffm1_fit$par)) {
ffm1_fit$par[decay_var] <- original_fit$par[decay_var]
} else {
ffm1_fit$fixed[decay_var,"value"] <- original_fit$fixed[decay_var,"value"]
}
new_col_name <- paste(metabolite, "DT50_FITTED_FFM_1", sep="_")
final[, new_col_name] <- CakeSinglePlot(ffm1_fit, xlim)[metabolite]
}
# Finally reset the fit for the next metabolite
fit <- original_fit
}
}
cat("<PLOT MODEL START>\n")
write.table(final, quote=FALSE)
cat("<PLOT MODEL END>\n")
# View(final)
}
# Get the immediate parents of a metabolite in a topology
GetParents <- function(topology, metabolite)
{
parents <- NULL
for (topo_metab in names(topology)) {
if (metabolite %in% topology[[topo_metab]]$to) {
parents <- append(parents, topo_metab)
}
}
return (parents)
}
# Get all the ancestors of a metabolite in a topology
GetAncestors <- function(topology, metabolite, ancestors = NULL)
{
for (parent in GetParents(topology, metabolite)) {
if (!(parent %in% ancestors)) {
ancestors <- GetAncestors(topology, parent, append(ancestors, parent))
}
}
return (ancestors)
}
# Set the formation leading to the metabolite to 1 from each parent
# Returns the updated fit or NULL if the metabolite can't be updated for FFM=1
SetFormationFractionsToOne <- function(fit, metabolite)
{
parents <- GetParents(fit$topology, metabolite)
for (parent in parents) {
if (any(GetAncestors(fit$topology, parent) %in% parents)) {
# Reject metabolites where any parent is also an ancestor of another parent
return (NULL)
}
# Try to find the formation fraction in the fitted and fixed parameters
ffm_var <- paste("f", parent, "to", metabolite, sep="_")
updated_fit <- SetFfmVariable(fit, ffm_var, 1)
if (!is.null(updated_fit)) {
# Formation fraction was found and set to 1
fit <- updated_fit
next
}
# The formation fraction for this metabolite is implicitly 1 minus all the others from this parent
# We need to iterate other formation fractions from the parent and set them to zero
for (other_metabolite in fit$topology[[parent]]$to) {
if (other_metabolite == metabolite) {
# Skip the metabolite itself because we know there isn't a formation fraction for it
next
}
other_ffm_var <- paste("f", parent, "to", other_metabolite, sep="_")
updated_fit <- SetFfmVariable(fit, other_ffm_var, 0)
# The updated fit shouldn't be able to be null because otherwise the topology was invalid
# or didn't match the parameters in the fit. If it does, we'll throw an error.
if (is.null(updated_fit)) {
stop(paste("The updated fit was NULL. The formation fraction '", other_ffm_var, "' doesn't exist in the fit.", sep=""))
}
fit <- updated_fit
}
}
return (fit)
}
# Set the value of ffm_var in the fit, looking in both fitted and fixed parameters
# Returns the updated fit or NULL if the ffm_var couldn't be found
SetFfmVariable <- function(fit, ffm_var, value)
{
if (ffm_var %in% names(fit$par)) {
# The ffm_var was found as a fitted parameter
fit$par[ffm_var] <- value
return (fit)
}
if (ffm_var %in% rownames(fit$fixed)) {
# The ffm_var was found as a fixed parameter
fit$fixed[ffm_var,"value"] <- value
return (fit)
}
# ffm_var was not found in the fit
return (NULL)
}
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-5
}
outtimes <- seq(0, xlim[2], length.out=CAKE.plots.resolution) * scale_x
# Solve the system
odeResult <- CakeOdeSolve(fit, outtimes, solution, atol)
odeResult_transformed <- PostProcessOdeOutput(odeResult, fit, atol)
# Replace values that are incalculably small with 0.
for (compartment.name in names(fit$map)) {
if (length(odeResult_transformed[, compartment.name][!is.nan(odeResult_transformed[, compartment.name])]) > 0) {
odeResult_transformed[, compartment.name][is.nan(odeResult_transformed[, compartment.name])] <- 0
}
odeResult_transformed[, compartment.name][odeResult_transformed[, compartment.name] < 0] <- 0
}
return (odeResult_transformed)
}
