Version 1.4v1.4

1 files changed, 117 insertions, 0 deletions

diff --git a/CakeOlsPlot.R b/CakeOlsPlot.R
new file mode 100644
index 0000000..199cc28
--- /dev/null
+++ b/CakeOlsPlot.R
@@ -0,0 +1,117 @@
+#$Id: CakeOlsPlot.R 216 2011-07-05 14:35:03Z nelr $
+# Generates fitted curves so the GUI can plot them
+# Based on code in IRLSkinfit
+# Author: Rob Nelson (Tessella)
+# Modifications developed by Tessella Plc for Syngenta: Copyright (C) 2011  Syngenta
+# Tessella Project Reference: 6245
+#
+#    This program is 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/>.”
+
+CakeOlsPlot <- function(fit, xlim = range(fit$data$time), ...)
+{
+   cat("<PLOT MODEL START>\n")
+
+   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) <- names(fit$diffs)
+
+  outtimes <- seq(0, xlim[2], length.out=101)
+
+  odenames <- c(
+    rownames(subset(fit$start, type == "deparm")),
+    rownames(subset(fit$fixed, type == "deparm")))
+  odeparms <- parms.all[odenames]
+
+  # 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("k1"), evalparse("k2"),
+          evalparse("g")),
+      HS = HS.solution(outtimes,
+          evalparse(parent.name),
+          evalparse("k1"), evalparse("k2"),
+          evalparse("tb")),
+      SFORB = SFORB.solution(outtimes,
+          evalparse(parent.name),
+          evalparse(paste("k", parent.name, "free_bound", sep="_")),
+          evalparse(paste("k", parent.name, "bound_free", sep="_")),
+          evalparse(paste("k", parent.name, "free_sink", sep="_")))
+    )
+    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
+    )
+  }
+    
+  # Output transformation for models with unobserved compartments like SFORB
+  out_transformed <- data.frame(time = out[,"time"])
+  for (var in names(fit$map)) {
+    if(length(fit$map[[var]]) == 1) {
+      out_transformed[var] <- out[, var]
+    } else {
+      out_transformed[var] <- rowSums(out[, fit$map[[var]]])
+    }
+  }
+  print(out_transformed)
+
+  cat("<PLOT MODEL END>\n")
+}