diff options
Diffstat (limited to 'CakeOlsPlot.R')
| -rw-r--r-- | CakeOlsPlot.R | 117 |
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") +} |