summaryrefslogtreecommitdiff
path: root/CakeCost.R
diff options
context:
space:
mode:
Diffstat (limited to 'CakeCost.R')
-rw-r--r--CakeCost.R225
1 files changed, 103 insertions, 122 deletions
diff --git a/CakeCost.R b/CakeCost.R
index 40635f0..8fb94ef 100644
--- a/CakeCost.R
+++ b/CakeCost.R
@@ -1,15 +1,14 @@
## -----------------------------------------------------------------------------
## The model cost and residuals
## -----------------------------------------------------------------------------
-# The CAKE R modules are based on mkin,
+# Some of the CAKE R modules are based on mkin.
# Call to approx is only performed if there are multiple non NA values
# which should prevent most of the crashes - Rob Nelson (Tessella)
#
-# Modifications developed by Tessella Plc for Syngenta, Copyright (C) 2011 Syngenta
-# Authors: Rob Nelson, Richard Smith
-# Tessella Project Reference: 6245
+# Modifications developed by Tessella for Syngenta, Copyright (C) 2011-2016 Syngenta
+# Tessella Project Reference: 6245, 7247, 8361, 7414
#
-# This program is free software: you can
+# The CAKE R modules are free software: you can
# redistribute them and/or modify them 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
@@ -21,9 +20,7 @@
# 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/>#
-#
-# $Id$
+# this program. If not, see <http://www.gnu.org/licenses/>
CakeCost <- function (model, obs, x = "time", y = NULL, err = NULL,
weight = "none", scaleVar = FALSE, cost = NULL, ...) {
@@ -180,23 +177,6 @@ CakeCost <- function (model, obs, x = "time", y = NULL, err = NULL,
Res <- (ModVar- obsdat)
res <- Res/Err
resScaled <- res*Scale
-
- # print("===========================");
- # print("#Names[i]");
- # print(Names[i]);
- # print("xDat");
- # print(xDat);
- # print("obsdat");
- # print(obsdat);
- # print("ModVar");
- # print(ModVar);
- # print("1/Err");
- # print(1/Err);
- # print("Res");
- # print(Res);
- # print("res");
- # print(res);
- # print("===========================");
Residual <- rbind(Residual,
data.frame(
@@ -217,9 +197,6 @@ CakeCost <- function (model, obs, x = "time", y = NULL, err = NULL,
SSR.unscaled = sum(res^2),
SSR = sum(resScaled^2)))
- # print("************************")
- # print(CostVar)
- # print("************************")
} # end loop over all observed variables
## SSR
@@ -268,103 +245,107 @@ plot.modCost<- function(x, legpos="topleft", ...) {
## Internal cost function for optimisers
## -----------------------------------------------------------------------------
# Cost function. The returned structure must have $model
-# Called from the middle of IrlsFit and OlsFit
# We need to preserve state between calls so make a closure
CakeInternalCostFunctions <- function(mkinmod, state.ini.optim, state.ini.optim.boxnames,
- state.ini.fixed, parms.fixed, observed, mkindiff, scaleVar,
- quiet, plot=FALSE, atol=1e-6){
+ state.ini.fixed, parms.fixed, observed, mkindiff,
+ quiet, atol=1e-6, solution="deSolve", err="err"){
cost.old <- 1e+100
calls <- 0
out_predicted <- NA
-
- get.predicted <- function(){ out_predicted }
-
- get.best.cost <- function(){ cost.old }
- reset.best.cost <- function() { cost.old<<-1e+100 }
-
- get.calls <- function(){ calls }
- set.calls <- function(newcalls){ calls <<- newcalls }
-
- set.error<-function(err) { observed$err <<- err }
-
- cost <- function(P) {
- assign("calls", calls + 1, inherits = TRUE)
- print(P)
- if (length(state.ini.optim) > 0) {
- odeini <- c(P[1:length(state.ini.optim)], state.ini.fixed)
- names(odeini) <- c(state.ini.optim.boxnames, names(state.ini.fixed))
- }
- else odeini <- state.ini.fixed
- odeparms <- c(P[(length(state.ini.optim) + 1):length(P)],
- parms.fixed)
- # Ensure initial state is at time 0
- outtimes = unique(c(0,observed$time))
- out <- ode(y = odeini, times = outtimes, func = mkindiff,
- parms = odeparms, atol=atol)
- out_transformed <- data.frame(time = out[, "time"])
- for (var in names(mkinmod$map)) {
- if (length(mkinmod$map[[var]]) == 1) {
- out_transformed[var] <- out[, var]
- }
- else {
- out_transformed[var] <- rowSums(out[, mkinmod$map[[var]]])
- }
- }
- assign("out_predicted", out_transformed, inherits = TRUE)
- mC <- CakeCost(out_transformed, observed, y = "value",
- err = 'err', scaleVar = scaleVar)
- mC$penalties <- CakePenalties(odeparms, out_transformed, observed)
- mC$model <- mC$cost + mC$penalties;
- if (mC$model < cost.old) {
- if (!quiet)
- cat("Model cost at call ", calls, ": m", mC$cost, 'p:', mC$penalties, 'o:', mC$model,
- "\n")
- if (plot) {
- outtimes_plot = seq(min(observed$time), max(observed$time),
- length.out = 100)
- out_plot <- ode(y = odeini, times = outtimes_plot,
- func = mkindiff, parms = odeparms, atol=atol)
- out_transformed_plot <- data.frame(time = out_plot[,
- "time"])
- for (var in names(mkinmod$map)) {
- if (length(mkinmod$map[[var]]) == 1) {
- out_transformed_plot[var] <- out_plot[, var]
- }
- else {
- out_transformed_plot[var] <- rowSums(out_plot[,
- mkinmod$map[[var]]])
- }
- }
- plot(0, type = "n", xlim = range(observed$time),
- ylim = range(observed$value, na.rm = TRUE),
- xlab = "Time", ylab = "Observed")
- col_obs <- pch_obs <- 1:length(obs_vars)
- names(col_obs) <- names(pch_obs) <- obs_vars
- for (obs_var in obs_vars) {
- points(subset(observed, name == obs_var, c(time,
- value)), pch = pch_obs[obs_var], col = col_obs[obs_var])
- }
- matlines(out_transformed_plot$time, out_transformed_plot[-1])
- legend("topright", inset = c(0.05, 0.05), legend = obs_vars,
- col = col_obs, pch = pch_obs, lty = 1:length(pch_obs))
- }
- assign("cost.old", mC$model, inherits = TRUE)
- }
- # HACK to make nls.lm respect the penalty, as it just uses residuals and ignores the cost
- if(mC$penalties > 0){
- #cat("Penalty adjustment: ",mC$penalties)
- mC$residuals$res <- mC$residuals$res + mC$penalties / length(mC$residuals$res)
- }
-
- return(mC)
- }
-
-
-
- list(cost=cost,
- get.predicted=get.predicted,
- get.calls=get.calls, set.calls=set.calls,
- get.best.cost=get.best.cost, reset.best.cost=reset.best.cost,
- set.error=set.error
- )
+
+ get.predicted <- function(){ out_predicted }
+
+ get.best.cost <- function(){ cost.old }
+ reset.best.cost <- function() { cost.old<<-1e+100 }
+
+ get.calls <- function(){ calls }
+ set.calls <- function(newcalls){ calls <<- newcalls }
+
+ set.error<-function(err) { observed$err <<- err }
+
+ # The called cost function
+ cost <- function(P) {
+ assign("calls", calls + 1, inherits = TRUE)
+ print(P)
+
+ if (length(state.ini.optim) > 0) {
+ odeini <- c(P[1:length(state.ini.optim)], state.ini.fixed)
+ names(odeini) <- c(state.ini.optim.boxnames, names(state.ini.fixed))
+ } else {
+ odeini <- state.ini.fixed
+ }
+
+ odeparms <- c(P[(length(state.ini.optim) + 1):length(P)], parms.fixed)
+
+ # Ensure initial state is at time 0
+ outtimes = unique(c(0,observed$time))
+
+ odeini <- AdjustOdeInitialValues(odeini, mkinmod, odeparms)
+
+ if (solution == "analytical") {
+ evalparse <- function(string)
+ {
+ eval(parse(text=string), as.list(c(odeparms, odeini)))
+ }
+
+ parent.type = names(mkinmod$map[[1]])[1]
+ parent.name = names(mkinmod$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", sub("_free", "", parent.name)))
+ }
+ if (solution == "deSolve")
+ {
+ out <- ode(y = odeini, times = outtimes, func = mkindiff, parms = odeparms, atol = atol)
+ }
+
+ out_transformed <- PostProcessOdeOutput(out, mkinmod, atol)
+
+ assign("out_predicted", out_transformed, inherits = TRUE)
+ mC <- CakeCost(out_transformed, observed, y = "value", err = err, scaleVar = FALSE)
+ mC$penalties <- CakePenalties(odeparms, out_transformed, observed)
+ mC$model <- mC$cost + mC$penalties
+
+ if (mC$model < cost.old) {
+ if (!quiet) {
+ cat("Model cost at call ", calls, ": m", mC$cost, 'p:', mC$penalties, 'o:', mC$model, "\n")
+ }
+
+ assign("cost.old", mC$model, inherits = TRUE)
+ }
+
+ # HACK to make nls.lm respect the penalty, as it just uses residuals and ignores the cost
+ if(mC$penalties > 0){
+ mC$residuals$res <- mC$residuals$res + mC$penalties / length(mC$residuals$res)
+ }
+
+ return(mC)
+ }
+
+ list(cost=cost,
+ get.predicted=get.predicted,
+ get.calls=get.calls, set.calls=set.calls,
+ get.best.cost=get.best.cost, reset.best.cost=reset.best.cost,
+ set.error=set.error
+ )
} \ No newline at end of file

Contact - Imprint