Version 3.2v3.2

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)
+}