diff options
Diffstat (limited to 'CakeCost.R')
| -rw-r--r-- | CakeCost.R | 225 |
1 files changed, 103 insertions, 122 deletions
@@ -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 |